<!DOCTYPE html>

<html>

<head>

<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />


<meta name="author" content="Ye Zheng, Kami Ahmad, Steven Henikoff" />


<title>CUT&amp;Tag Data Processing and Analysis Tutorial</title>

<script src="site_libs/jquery-1.11.3/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/cosmo.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="site_libs/jqueryui-1.11.4/jquery-ui.min.js"></script>
<link href="site_libs/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="site_libs/tocify-1.9.1/jquery.tocify.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<script src="site_libs/navigation-1.1/codefolding.js"></script>
<link href="site_libs/highlightjs-9.12.0/textmate.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<link href="site_libs/pagedtable-1.1/css/pagedtable.css" rel="stylesheet" />
<script src="site_libs/pagedtable-1.1/js/pagedtable.js"></script>
<link href="site_libs/font-awesome-5.1.0/css/all.css" rel="stylesheet" />
<link href="site_libs/font-awesome-5.1.0/css/v4-shims.css" rel="stylesheet" />

<link rel="icon" href="https://github.com/workflowr/workflowr-assets/raw/master/img/reproducible.png">
<!-- Add a small amount of space between sections. -->
<style type="text/css">
div.section {
  padding-top: 12px;
}
</style>

<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
  pre:not([class]) {
    background-color: white;
  }
</style>
<script type="text/javascript">
if (window.hljs) {
  hljs.configure({languages: []});
  hljs.initHighlightingOnLoad();
  if (document.readyState && document.readyState === "complete") {
    window.setTimeout(function() { hljs.initHighlighting(); }, 0);
  }
}
</script>



<style type="text/css">
h1 {
  font-size: 34px;
}
h1.title {
  font-size: 38px;
}
h2 {
  font-size: 30px;
}
h3 {
  font-size: 24px;
}
h4 {
  font-size: 18px;
}
h5 {
  font-size: 16px;
}
h6 {
  font-size: 12px;
}
.table th:not([align]) {
  text-align: left;
}
</style>




<style type = "text/css">
.main-container {
  max-width: 940px;
  margin-left: auto;
  margin-right: auto;
}
code {
  color: inherit;
  background-color: rgba(0, 0, 0, 0.04);
}
img {
  max-width:100%;
}
.tabbed-pane {
  padding-top: 12px;
}
.html-widget {
  margin-bottom: 20px;
}
button.code-folding-btn:focus {
  outline: none;
}
summary {
  display: list-item;
}
</style>


<style type="text/css">
/* padding for bootstrap navbar */
body {
  padding-top: 51px;
  padding-bottom: 40px;
}
/* offset scroll position for anchor links (for fixed navbar)  */
.section h1 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h2 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h3 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h4 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h5 {
  padding-top: 56px;
  margin-top: -56px;
}
.section h6 {
  padding-top: 56px;
  margin-top: -56px;
}
.dropdown-submenu {
  position: relative;
}
.dropdown-submenu>.dropdown-menu {
  top: 0;
  left: 100%;
  margin-top: -6px;
  margin-left: -1px;
  border-radius: 0 6px 6px 6px;
}
.dropdown-submenu:hover>.dropdown-menu {
  display: block;
}
.dropdown-submenu>a:after {
  display: block;
  content: " ";
  float: right;
  width: 0;
  height: 0;
  border-color: transparent;
  border-style: solid;
  border-width: 5px 0 5px 5px;
  border-left-color: #cccccc;
  margin-top: 5px;
  margin-right: -10px;
}
.dropdown-submenu:hover>a:after {
  border-left-color: #ffffff;
}
.dropdown-submenu.pull-left {
  float: none;
}
.dropdown-submenu.pull-left>.dropdown-menu {
  left: -100%;
  margin-left: 10px;
  border-radius: 6px 0 6px 6px;
}
</style>

<script>
// manage active state of menu based on current page
$(document).ready(function () {
  // active menu anchor
  href = window.location.pathname
  href = href.substr(href.lastIndexOf('/') + 1)
  if (href === "")
    href = "index.html";
  var menuAnchor = $('a[href="' + href + '"]');

  // mark it active
  menuAnchor.parent().addClass('active');

  // if it's got a parent navbar menu mark it active as well
  menuAnchor.closest('li.dropdown').addClass('active');
});
</script>

<!-- tabsets -->

<style type="text/css">
.tabset-dropdown > .nav-tabs {
  display: inline-table;
  max-height: 500px;
  min-height: 44px;
  overflow-y: auto;
  background: white;
  border: 1px solid #ddd;
  border-radius: 4px;
}

.tabset-dropdown > .nav-tabs > li.active:before {
  content: "";
  font-family: 'Glyphicons Halflings';
  display: inline-block;
  padding: 10px;
  border-right: 1px solid #ddd;
}

.tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before {
  content: "&#xe258;";
  border: none;
}

.tabset-dropdown > .nav-tabs.nav-tabs-open:before {
  content: "";
  font-family: 'Glyphicons Halflings';
  display: inline-block;
  padding: 10px;
  border-right: 1px solid #ddd;
}

.tabset-dropdown > .nav-tabs > li.active {
  display: block;
}

.tabset-dropdown > .nav-tabs > li > a,
.tabset-dropdown > .nav-tabs > li > a:focus,
.tabset-dropdown > .nav-tabs > li > a:hover {
  border: none;
  display: inline-block;
  border-radius: 4px;
  background-color: transparent;
}

.tabset-dropdown > .nav-tabs.nav-tabs-open > li {
  display: block;
  float: none;
}

.tabset-dropdown > .nav-tabs > li {
  display: none;
}
</style>

<!-- code folding -->
<style type="text/css">
.code-folding-btn { margin-bottom: 4px; }
</style>



<style type="text/css">

#TOC {
  margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
  position: relative;
  width: 100%;
}
}

@media print {
.toc-content {
  /* see https://github.com/w3c/csswg-drafts/issues/4434 */
  float: right;
}
}

.toc-content {
  padding-left: 30px;
  padding-right: 40px;
}

div.main-container {
  max-width: 1200px;
}

div.tocify {
  width: 20%;
  max-width: 260px;
  max-height: 85%;
}

@media (min-width: 768px) and (max-width: 991px) {
  div.tocify {
    width: 25%;
  }
}

@media (max-width: 767px) {
  div.tocify {
    width: 100%;
    max-width: none;
  }
}

.tocify ul, .tocify li {
  line-height: 20px;
}

.tocify-subheader .tocify-item {
  font-size: 0.90em;
}

.tocify .list-group-item {
  border-radius: 0px;
}


</style>



</head>

<body>


<div class="container-fluid main-container">


<!-- setup 3col/9col grid for toc_float and main content  -->
<div class="row-fluid">
<div class="col-xs-12 col-sm-4 col-md-3">
<div id="TOC" class="tocify">
</div>
</div>

<div class="toc-content col-xs-12 col-sm-8 col-md-9">




<div class="navbar navbar-default  navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <a class="navbar-brand" href="index.html">CUTTag_tutorial</a>
    </div>
    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
        <li>
  <a href="index.html">Home</a>
</li>
<li>
  <a href="about.html">About</a>
</li>
<li>
  <a href="license.html">License</a>
</li>
      </ul>
      <ul class="nav navbar-nav navbar-right">
        <li>
  <a href="https://github.com/yezhengSTAT/CUTTag_tutorial">
    <span class="fa fa-github"></span>
     
    Source code
  </a>
</li>
      </ul>
    </div><!--/.nav-collapse -->
  </div><!--/.container -->
</div><!--/.navbar -->

<div class="fluid-row" id="header">

<div class="btn-group pull-right">
<button type="button" class="btn btn-default btn-xs dropdown-toggle" data-toggle="dropdown" aria-haspopup="true" aria-expanded="false"><span>Code</span> <span class="caret"></span></button>
<ul class="dropdown-menu" style="min-width: 50px;">
<li><a id="rmd-show-all-code" href="#">Show All Code</a></li>
<li><a id="rmd-hide-all-code" href="#">Hide All Code</a></li>
</ul>
</div>



<h1 class="title toc-ignore">CUT&amp;Tag Data Processing and Analysis Tutorial</h1>
<h4 class="author">Ye Zheng, Kami Ahmad, Steven Henikoff</h4>

</div>


<p>
<button type="button" class="btn btn-default btn-workflowr btn-workflowr-report" data-toggle="collapse" data-target="#workflowr-report">
<span class="glyphicon glyphicon-list" aria-hidden="true"></span> workflowr <span class="glyphicon glyphicon-exclamation-sign text-danger" aria-hidden="true"></span>
</button>
</p>
<div id="workflowr-report" class="collapse">
<ul class="nav nav-tabs">
<li class="active">
<a data-toggle="tab" href="#summary">Summary</a>
</li>
<li>
<a data-toggle="tab" href="#checks"> Checks <span class="glyphicon glyphicon-exclamation-sign text-danger" aria-hidden="true"></span> </a>
</li>
<li>
<a data-toggle="tab" href="#versions">Past versions</a>
</li>
</ul>
<div class="tab-content">
<div id="summary" class="tab-pane fade in active">
<p>
<strong>Last updated:</strong> 2021-02-01
</p>
<p>
<strong>Checks:</strong> <span class="glyphicon glyphicon-ok text-success" aria-hidden="true"></span> 6 <span class="glyphicon glyphicon-exclamation-sign text-danger" aria-hidden="true"></span> 1
</p>
<p>
<strong>Knit directory:</strong> <code>CUTTag_tutorial/</code> <span class="glyphicon glyphicon-question-sign" aria-hidden="true" title="This is the local directory in which the code in this file was executed."> </span>
</p>
<p>
This reproducible <a href="http://rmarkdown.rstudio.com">R Markdown</a> analysis was created with <a
  href="https://github.com/jdblischak/workflowr">workflowr</a> (version 1.6.2). The <em>Checks</em> tab describes the reproducibility checks that were applied when the results were created. The <em>Past versions</em> tab lists the development history.
</p>
<hr>
</div>
<div id="checks" class="tab-pane fade">
<div class="panel-group" id="workflowr-checks">
<div class="panel panel-default">
<div class="panel-heading">
<p class="panel-title">
<a data-toggle="collapse" data-parent="#workflowr-checks" href="#strongRMarkdownfilestronguptodate"> <span class="glyphicon glyphicon-ok text-success" aria-hidden="true"></span> <strong>R Markdown file:</strong> up-to-date </a>
</p>
</div>
<div id="strongRMarkdownfilestronguptodate" class="panel-collapse collapse">
<div class="panel-body">
<p>Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.</p>
</div>
</div>
</div>
<div class="panel panel-default">
<div class="panel-heading">
<p class="panel-title">
<a data-toggle="collapse" data-parent="#workflowr-checks" href="#strongEnvironmentstrongempty"> <span class="glyphicon glyphicon-ok text-success" aria-hidden="true"></span> <strong>Environment:</strong> empty </a>
</p>
</div>
<div id="strongEnvironmentstrongempty" class="panel-collapse collapse">
<div class="panel-body">
<p>Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.</p>
</div>
</div>
</div>
<div class="panel panel-default">
<div class="panel-heading">
<p class="panel-title">
<a data-toggle="collapse" data-parent="#workflowr-checks" href="#strongSeedstrongcodesetseed20200415code"> <span class="glyphicon glyphicon-ok text-success" aria-hidden="true"></span> <strong>Seed:</strong> <code>set.seed(20200415)</code> </a>
</p>
</div>
<div id="strongSeedstrongcodesetseed20200415code" class="panel-collapse collapse">
<div class="panel-body">
<p>The command <code>set.seed(20200415)</code> was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.</p>
</div>
</div>
</div>
<div class="panel panel-default">
<div class="panel-heading">
<p class="panel-title">
<a data-toggle="collapse" data-parent="#workflowr-checks" href="#strongSessioninformationstrongrecorded"> <span class="glyphicon glyphicon-ok text-success" aria-hidden="true"></span> <strong>Session information:</strong> recorded </a>
</p>
</div>
<div id="strongSessioninformationstrongrecorded" class="panel-collapse collapse">
<div class="panel-body">
<p>Great job! Recording the operating system, R version, and package versions is critical for reproducibility.</p>
</div>
</div>
</div>
<div class="panel panel-default">
<div class="panel-heading">
<p class="panel-title">
<a data-toggle="collapse" data-parent="#workflowr-checks" href="#strongCachestrongnone"> <span class="glyphicon glyphicon-ok text-success" aria-hidden="true"></span> <strong>Cache:</strong> none </a>
</p>
</div>
<div id="strongCachestrongnone" class="panel-collapse collapse">
<div class="panel-body">
<p>Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.</p>
</div>
</div>
</div>
<div class="panel panel-default">
<div class="panel-heading">
<p class="panel-title">
<a data-toggle="collapse" data-parent="#workflowr-checks" href="#strongFilepathsstrongabsolute"> <span class="glyphicon glyphicon-exclamation-sign text-danger" aria-hidden="true"></span> <strong>File paths:</strong> absolute </a>
</p>
</div>
<div id="strongFilepathsstrongabsolute" class="panel-collapse collapse">
<div class="panel-body">
<p>
Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.
</p>
<table class="table table-condensed table-hover">
<thead>
<tr>
<th style="text-align:left;">
absolute
</th>
<th style="text-align:left;">
relative
</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align:left;">
/fh/fast/gottardo_r/yezheng_working/cuttag/CUTTag_tutorial
</td>
<td style="text-align:left;">
.
</td>
</tr>
</tbody>
</table>
</div>
</div>
</div>
<div class="panel panel-default">
<div class="panel-heading">
<p class="panel-title">
<a data-toggle="collapse" data-parent="#workflowr-checks" href="#strongRepositoryversionstrongahrefhttpsgithubcomyezhengSTATCUTTagtutorialtree4ff93d77d1a1dfd05ff69489d428c098a1779703targetblank4ff93d7a"> <span class="glyphicon glyphicon-ok text-success" aria-hidden="true"></span> <strong>Repository version:</strong> <a href="https://github.com/yezhengSTAT/CUTTag_tutorial/tree/4ff93d77d1a1dfd05ff69489d428c098a1779703" target="_blank">4ff93d7</a> </a>
</p>
</div>
<div id="strongRepositoryversionstrongahrefhttpsgithubcomyezhengSTATCUTTagtutorialtree4ff93d77d1a1dfd05ff69489d428c098a1779703targetblank4ff93d7a" class="panel-collapse collapse">
<div class="panel-body">
<p>
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
</p>
<p>
The results in this page were generated with repository version <a href="https://github.com/yezhengSTAT/CUTTag_tutorial/tree/4ff93d77d1a1dfd05ff69489d428c098a1779703" target="_blank">4ff93d7</a>. See the <em>Past versions</em> tab to see a history of the changes made to the R Markdown and HTML files.
</p>
<p>
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use <code>wflow_publish</code> or <code>wflow_git_commit</code>). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
</p>
<pre><code>
Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/.DS_Store

Untracked files:
    Untracked:  ._.DS_Store
    Untracked:  alignment/
    Untracked:  data/IgG_old/
    Untracked:  data/IgG_rep1/
    Untracked:  data/IgG_rep2/
    Untracked:  data/K27me3_rep1/
    Untracked:  data/K27me3_rep2/
    Untracked:  data/K4me3_rep1/
    Untracked:  data/K4me3_rep2/
    Untracked:  data/hg38_gene/
    Untracked:  fastq/
    Untracked:  fastqFileQC/
    Untracked:  peakCalling/

Unstaged changes:
    Modified:   analysis/pipeline2.sh

</code></pre>
<p>
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
</p>
</div>
</div>
</div>
</div>
<hr>
</div>
<div id="versions" class="tab-pane fade">

<p>
These are the previous versions of the repository in which changes were made to the R Markdown (<code>analysis/tutorials.Rmd</code>) and HTML (<code>docs/tutorials.html</code>) files. If you’ve configured a remote Git repository (see <code>?wflow_git_remote</code>), click on the hyperlinks in the table below to view the files as they were in that past version.
</p>
<div class="table-responsive">
<table class="table table-condensed table-hover">
<thead>
<tr>
<th>
File
</th>
<th>
Version
</th>
<th>
Author
</th>
<th>
Date
</th>
<th>
Message
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
Rmd
</td>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/4ff93d77d1a1dfd05ff69489d428c098a1779703/analysis/tutorials.Rmd" target="_blank">4ff93d7</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2021-02-01
</td>
<td>
edit typo
</td>
</tr>
<tr>
<td>
Rmd
</td>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/421e40f4713aff3af285e6939b1f55dd407633e6/analysis/tutorials.Rmd" target="_blank">421e40f</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-20
</td>
<td>
update minor issues
</td>
</tr>
<tr>
<td>
html
</td>
<td>
<a href="https://rawcdn.githack.com/yezhengSTAT/CUTTag_tutorial/b2993411ac609462d4ebc4b7ce1446e327a00e56/docs/tutorials.html" target="_blank">b299341</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
<td>
Build site.
</td>
</tr>
<tr>
<td>
Rmd
</td>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/438c1a5881735ee1c4e4ebe25542173b07a9cf98/analysis/tutorials.Rmd" target="_blank">438c1a5</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
<td>
Publish the initial files for myproject
</td>
</tr>
<tr>
<td>
Rmd
</td>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/d7f438eb1c959c4b627814b1aa1f07850756b6bd/analysis/tutorials.Rmd" target="_blank">d7f438e</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
<td>
update minor issues
</td>
</tr>
<tr>
<td>
html
</td>
<td>
<a href="https://rawcdn.githack.com/yezhengSTAT/CUTTag_tutorial/d7f438eb1c959c4b627814b1aa1f07850756b6bd/docs/tutorials.html" target="_blank">d7f438e</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
<td>
update minor issues
</td>
</tr>
<tr>
<td>
Rmd
</td>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/1118ace8ae76a1dfa0d38db6fd0b79767d0d66ec/analysis/tutorials.Rmd" target="_blank">1118ace</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-07-16
</td>
<td>
update tutorial
</td>
</tr>
<tr>
<td>
Rmd
</td>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/d443a2492f74c43e61305bc13f5e0e54183388e8/analysis/tutorials.Rmd" target="_blank">d443a24</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-06-01
</td>
<td>
add figure
</td>
</tr>
<tr>
<td>
html
</td>
<td>
<a href="https://rawcdn.githack.com/yezhengSTAT/CUTTag_tutorial/8db752dee6286a2e0931b1ba6832bc06868d37b7/docs/tutorials.html" target="_blank">8db752d</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-06-01
</td>
<td>
Build site.
</td>
</tr>
<tr>
<td>
html
</td>
<td>
<a href="https://rawcdn.githack.com/yezhengSTAT/CUTTag_tutorial/2e55365abe56f321e185f5bf0262114c9484e658/docs/tutorials.html" target="_blank">2e55365</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-06-01
</td>
<td>
Build site.
</td>
</tr>
<tr>
<td>
html
</td>
<td>
<a href="https://rawcdn.githack.com/yezhengSTAT/CUTTag_tutorial/d544756d9ed67d5e7ab554eb1fd2fc9e410c5f13/docs/tutorials.html" target="_blank">d544756</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-06-01
</td>
<td>
Build site.
</td>
</tr>
<tr>
<td>
Rmd
</td>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/77afcb1c79ef7b513e7ea9617e5837118552de3b/analysis/tutorials.Rmd" target="_blank">77afcb1</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-06-01
</td>
<td>
Publish the initial files for myproject
</td>
</tr>
</tbody>
</table>
</div>
<hr>
</div>
</div>
</div>
<!-- --- -->
<!-- title: "CUT&Tag Data Processing and Analysis Tutorial" -->
<!-- author: "Ye Zheng" -->
<!-- output: -->
<!--   html_document: -->
<!--     toc: true -->
<!--     toc_float: true -->
<!--     #number_sections: false -->
<!--     code_folding: show -->
<!-- editor_options: -->
<!--   chunk_output_type: console -->
<!-- --- -->
<!-- <style> -->
<!-- body{text-align: justify} -->
<!-- pre code, pre, code { -->
<!--   white-space: pre !important; -->
<!--   overflow-x: scroll !important; -->
<!--   word-break: keep-all !important; -->
<!--   word-wrap: initial !important; -->
<!-- } -->
<!-- </style> -->
<p><em>contact: <a href="mailto:yzheng23@fredhutch.org" class="email">yzheng23@fredhutch.org</a></em></p>
<p><em>File creation: March 13, 2020</em></p>
<p><em>Update: Aug 10, 2020</em></p>
<p><strong>Approximate time: 60 - 120 minutes</strong></p>
<div id="i.-introduction" class="section level1">
<h1>I. Introduction</h1>
<div id="overview-of-cuttag" class="section level2">
<h2>1.1. Overview of CUT&amp;Tag</h2>
<p>All dynamic processes that take place on DNA in the eukaryotic nucleus occur in the context of a chromatin landscape that comprises nucleosomes and their modifications, transcription factors, and chromatin-associated complexes. A variety of chromatin features mark sites of activating and silencing transcriptional regulatory elements and chromatin domains that differ between cell types and change during development.</p>
<p>The mapping of chromatin features genome-wide has traditionally been performed using chromatin immunoprecipitation (ChIP), in which chromatin is cross-linked and solubilized, and an antibody to a protein or modification of interest is used to immunoprecipitate the bound DNA (Fig. 1a). Very little has changed in the way ChIP is most generally performed since it was first described 35 years ago, and remains fraught with signal-to-noise issues and artifacts. An alternative chromatin profiling strategy is enzyme tethering in situ whereby the chromatin protein or modification of interest is targeted by an antibody or fusion protein. Then, the underlying DNA is marked or cleaved, and a succession of enzyme-tethering methods have been introduced over the past two decades. Cleavage Under Targets &amp; Tagmentation (CUT&amp;Tag) is a tethering method that uses a protein-A-Tn5 (pA-Tn5) transposome fusion protein (Fig. 1b). In CUT&amp;Tag, permeabilized cells or nuclei are incubated with antibody to a specified chromatin protein, and then pA-Tn5 loaded with mosaic end adaptors is successively tethered to antibody-bound sites. Activation of the transposome by adding magnesium ions results in the integration of the adaptors into nearby DNA. These are then amplified to generate sequencing libraries. Antibody-tethered Tn5-based methods achieve high sensitivity owing to stringent washing of samples after pA-Tn5 tethering and the high efficiency of adaptor integration. The improved signal-to-noise relative to ChIP-seq translates to an order-of-magnitude reduction in the amount of sequencing required to map chromatin features, allowing sample pooling (typically up to 90 samples) for paired-end sequencing on Illumina NGS sequencers by barcoded PCR of libraries.</p>
<div class="figure">
<img src="figures/ChIPseqCUTTag.png" alt="Figure 1. Differences between immunoprecipitation and in antibody-targeted chromatin profiling strategies. A. ChIP-seq experimental procedure. B. CUT&amp;Tag experimental procedure. Cells and nuclei are indicated in grey, chromatin as red nucleosomes, and a specific chromatin protein in green." />
<p class="caption"><strong>Figure 1. Differences between immunoprecipitation and in antibody-targeted chromatin profiling strategies.</strong> <strong>A.</strong> ChIP-seq experimental procedure. <strong>B.</strong> CUT&amp;Tag experimental procedure. Cells and nuclei are indicated in grey, chromatin as red nucleosomes, and a specific chromatin protein in green.</p>
</div>
</div>
<div id="objectives" class="section level2">
<h2>1.2. Objectives</h2>
<p>This tutorial is designed for processing and analyzing CUT&amp;Tag data following the <a href="https://www.protocols.io/view/bench-top-cut-amp-tag-bcuhiwt6/abstract">Bench top CUT&amp;Tag V.3 protocol</a>. The illustration data used in this tutorial is the profiling of histone modifications in the human lymphoma K562 cell line, but the tutorial is generally applicable to any chromatin protein, including transcription factors, RNA polymerase II, and epitope-tagged proteins.</p>
</div>
<div id="cuttag-data-processing-and-analysis-outline." class="section level2">
<h2>1.3 CUT&amp;Tag data processing and analysis outline.</h2>
<div class="figure">
<img src="figures/CUTTag_Tuto_Diagram.png" alt="Figure 2. CUT&amp;Tag data processing and analysis." />
<p class="caption"><strong>Figure 2. CUT&amp;Tag data processing and analysis.</strong></p>
</div>
</div>
<div id="requirements" class="section level2">
<h2>1.4. Requirements</h2>
<ul>
<li><p>Linux system</p></li>
<li>R (versions &gt;= 3.6)
<ul>
<li>dplyr</li>
<li>stringr</li>
<li>ggplot2</li>
<li>viridis</li>
<li>GenomicRanges</li>
<li>chromVAR</li>
<li>DESeq2</li>
<li>ggpubr</li>
<li>corrplot</li>
<li>ChIPseqSpikeInFree [Optional]</li>
</ul></li>
</ul>
<pre class="r"><code>library(dplyr)
library(stringr)
library(ggplot2)
library(viridis)
library(GenomicRanges)
library(chromVAR) ## For FRiP analysis and differential analysis
library(DESeq2) ## For differential analysis section
library(ggpubr) ## For customizing figures
library(corrplot) ## For correlation plot</code></pre>
<ul>
<li><p>FastQC(version &gt;= 0.11.9)</p></li>
<li><p>Bowtie2 (version &gt;= 2.3.4.3)</p></li>
<li><p>samtools (version &gt;= 1.10)</p></li>
<li><p>bedtools (version &gt;= 2.29.1)</p></li>
<li><p>Picard (version &gt;= 2.18.29)</p></li>
<li><p>SEACR (version &gt;= 1.3)</p></li>
<li><p>deepTools (version &gt;= 2.0)</p></li>
</ul>
</div>
<div id="data-downloading" class="section level2">
<h2>1.5. Data Downloading</h2>
<p>In this tutorial, we use data from Kaya-Okur et al. (2020), and available for downloading from <a href="https://www.ncbi.nlm.nih.gov/geo/">GEO</a>. The corresponding SRA entries are provided below.</p>
<p>Options to download SRA sequences from GEO:</p>
<ul>
<li><ol style="list-style-type: lower-alpha">
<li>Using <a href="https://www.ncbi.nlm.nih.gov/sra/docs/sradownload/">SRA Toolkit</a></li>
</ol></li>
<li><ol start="2" style="list-style-type: lower-alpha">
<li>Download through <a href="file:///Users/yezheng/Downloads/downloading_fastq_GEO.pdf">European Nucleotide Archive</a>. New ENA Browser: <a href="https://www.ebi.ac.uk/ena/browser/view" class="uri">https://www.ebi.ac.uk/ena/browser/view</a>. We are using this option as illustration.</li>
</ol></li>
</ul>
<p>First, we need to specify the project path.</p>
<pre class="bash"><code>##== linux command ==##
projPath=&quot;/path/to/project/where/data/and/results/are/saved&quot;</code></pre>
<ul>
<li>H3K27me3:
<ul>
<li>SH_Hs_K27m3_NX_0918 as replicate 1: GEO accession: GSE145187, SRA entry: SRX8754646</li>
<li>SH_Hs_K27m3_Xpc_0107 as replicate 2: GEO accession: GSE145187, SRA entry: SRX7713678</li>
</ul></li>
<li>H3K4me3:
<ul>
<li>SH_Hs_K4m3_NX_0918 as replicate 1: GEO accession: GSE145187, SRA entry: SRX7713692</li>
<li>SH_Hs_K4m3_Xpc_0107 as replicate 2: GEO accession: GSE145187, SRA entry: SRX7713696</li>
</ul></li>
<li>IgG:
<ul>
<li>SH_Hs_IgG_1x_0924 as replicate 1:GEO accession: GSE145187, SRA entry: SRX8468909</li>
<li>SH_Hs_IgG_20181224 as replicate 2: GEO accession: GSM3680227, SRA entry: SRX5545346</li>
</ul></li>
</ul>
<p>Taking SH_Hs_IgG_20181224 as an example.</p>
<pre class="bash"><code>##== linux command ==##
wget -O $projPath/data/IgG_rep2/IgG_rep2_R1_001.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/001/SRR8754611/SRR8754611_1.fastq.gz

wget -O $projPath/data/IgG_rep2/IgG_rep2_R2_001.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/001/SRR8754611/SRR8754611_2.fastq.gz

wget -O $projPath/data/IgG_rep2/IgG_rep2_R1_002.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/002/SRR8754612/SRR8754612_1.fastq.gz

wget -O $projPath/data/IgG_rep2/IgG_rep2_R2_002.fastq.gz ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR875/002/SRR8754612/SRR8754612_2.fastq.gz</code></pre>
</div>
</div>
<div id="ii.-data-pre-processing" class="section level1">
<h1>II. Data Pre-processing</h1>
<div id="quality-control-using-fastqc-optional" class="section level2">
<h2>2.1. Quality Control using <a href="https://www.bioinformatics.babraham.ac.uk/projects/fastqc/">FastQC</a> [Optional]</h2>
<p>This step is not required. In case that users are generating their own data and FastQC is one of the routine checking procedures in users’ group, we provide this step as troubleshoting explanation.</p>
<div id="obtain-fastqc" class="section level3">
<h3>2.2.1 Obtain FastQC</h3>
<pre class="bash"><code>##== linux command ==##
mkdir -p $projPath/tools
wget -P $projPath/tools https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.9.zip
cd $projPath/tools
unzip fastqc_v0.11.9.zip</code></pre>
</div>
<div id="run-fastqc-for-quality-check" class="section level3">
<h3>2.2.2 Run FastQC for quality check</h3>
<pre class="bash"><code>##== linux command ==##
mkdir -p ${projPath}/fastqFileQC/${histName}

$projPath/tools/FastQC/fastqc -o ${projPath}/fastqFileQC/${histName} -f fastq ${projPath}/fastq/${histName}_R1.fastq.gz
$projPath/tools/FastQC/fastqc -o ${projPath}/fastqFileQC/${histName} -f fastq ${projPath}/fastq/${histName}_R2.fastq.gz</code></pre>
</div>
<div id="intepret-the-quality-check-results." class="section level3">
<h3>2.2.3 Intepret the quality check results.</h3>
<p>Quality check reference: <a href="https://www.bioinformatics.babraham.ac.uk/projects/fastqc/bad_sequence_fastqc.html">https://www.bioinformatics.babraham.ac.uk/projects/fastqc/bad_sequence_fastqc.html</a></p>
<div class="figure">
<img src="figures/R1_sequenceContent.png" alt="Figure 3. Per base sequence content fails the FastQC quality check." />
<p class="caption"><strong>Figure 3. Per base sequence content fails the FastQC quality check.</strong></p>
</div>
<p><strong>The discordant sequence content at the begining of the reads are common phenomenon for CUT&amp;Tag reads. Failing to pass the Per base seuqnence content does not mean your data failed.</strong></p>
<ul>
<li><p>It can be due to the Tn5 preference.</p></li>
<li><p>What you might be detecting is the 10-bp periodicity that shows up as a sawtooth pattern in the length distribution. If so, this is normal and will not affect alignment or peak calling. In any case we do not recommend trimming as the bowtie2 parameters that we list will give accurate mapping information without trimming.</p></li>
</ul>
</div>
</div>
<div id="merge-technical-replicateslanes-if-needed-optional" class="section level2">
<h2>2.2. Merge technical replicates/lanes if needed [Optional]</h2>
<p>Sometimes, samples are often sequenced across multiple lanes for efficiency and can be pooled before alignment. If you want to check the reproducibility between sequences of different lanes of the same sample, you can skip this step and align each sequencing file (fastq file) respectively.</p>
<pre class="bash"><code>##== linux command ==##
histName=&quot;K27me3_rep1&quot;

mkdir -p ${projPath}/fastq
cat ${projPath}/data/${histName}/*_R1_*.fastq.gz &gt;${projPath}/fastq/${histName}_R1.fastq.gz
cat ${projPath}/data/${histName}/*_R2_*.fastq.gz &gt;${projPath}/fastq/${histName}_R2.fastq.gz</code></pre>
</div>
</div>
<div id="iii.-alignment" class="section level1">
<h1>III. Alignment</h1>
<div id="bowtie2-alignment-required" class="section level2">
<h2>3.1. Bowtie2 alignment [required]</h2>
<p>The structure of CUT&amp;Tag insert libraries with Tn5 adapters and barcoded PCR primers is shown below:</p>
<div class="figure">
<img src="figures/BarcodedCnTLibrary.png" alt="Figure 4. CUT&amp;Tag insert libraries with the sequence of adapters." />
<p class="caption"><strong>Figure 4. CUT&amp;Tag insert libraries with the sequence of adapters.</strong></p>
</div>
<p>Our standard pipeline is to perform single-index 25x25 PE Illumina sequencing on up to 90 pooled samples on a single HiSeq 2500 flowcell, where each sample has a unique PCR primer barcode. Amounts for each library are adjusted to provide ~5 million paired-end reads, which provides high-quality profiling for abundant chromatin features with a specific and high-yield antibody. Less abundant features typically require fewer reads, while lower-quality antibodies may increase the number of reads needed for generating robust chromatin profiles. A thorough discussion of feature recall and sequencing depths for CUT&amp;Tag has been published (Kaya-Okur et al 2020).</p>
<div id="alignment-to-hg38." class="section level3">
<h3>3.1.1 Alignment to HG38.</h3>
<pre class="bash"><code>##== linux command ==##
cores=8
ref=&quot;/path/to/bowtie2Index/hg38&quot;

mkdir -p ${projPath}/alignment/sam/bowtie2_summary
mkdir -p ${projPath}/alignment/bam
mkdir -p ${projPath}/alignment/bed
mkdir -p ${projPath}/alignment/bedgraph

## Build the bowtie2 reference genome index if needed:
## bowtie2-build path/to/hg38/fasta/hg38.fa /path/to/bowtie2Index/hg38

bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 -p ${cores} -x ${ref} -1 ${projPath}/fastq/${histName}_R1.fastq.gz -2 ${projPath}/fastq/${histName}_R2.fastq.gz -S ${projPath}/alignment/sam/${histName}_bowtie2.sam &amp;&gt; ${projPath}/alignment/sam/bowtie2_summary/${histName}_bowtie2.txt</code></pre>
<p>The paired-end reads are aligned by Bowtie2 using parameters <code>--end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700</code> for mapping of inserts 10-700 bp in length.</p>
<p><strong>Critical step</strong>: There is no need to trim reads from out standard 25x25 PE sequencing, as adapter sequences will not be included in reads of inserts &gt;25 bp. However, for users performing longer sequencing, reads will need to be trimmed by Cutadapt and mapped by <code>--local --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700</code> to ignore any remaining adapter sequence at the 3’ ends of reads during mapping.</p>
</div>
<div id="alignment-to-spike-in-genome-for-spike-in-calibration-optionalrecommended" class="section level3">
<h3>3.1.2 Alignment to spike-in genome for spike-in calibration [optional/recommended]</h3>
<p>This section is <strong>optional</strong> but <strong>recommended</strong> depending on your experimental protocol.</p>
<p>E. coli DNA is carried along with bacterially-produced pA-Tn5 protein and gets tagmented non-specifically during the reaction. The fraction of total reads that map to the E.coli genome depends on the yield of epitope-targeted CUT&amp;Tag, and so depends on the number of cells used and the abundance of that epitope in chromatin. Since a constant amount of pATn5 is added to CUT&amp;Tag reactions and brings along a fixed amount of E. coli DNA, E. coli reads can be used to normalize epitope abundance in a set of experiments. More discussion, please see Section V.</p>
<pre class="bash"><code>##== linux command ==##
spikeInRef=&quot;/shared/ngs/illumina/henikoff/Bowtie2/Ecoli&quot;
chromSize=&quot;/fh/fast/gottardo_r/yezheng_working/SupplementaryData/hg38/chromSize/hg38.chrom.size&quot;

## bowtie2-build path/to/Ecoli/fasta/Ecoli.fa /path/to/bowtie2Index/Ecoli
bowtie2 --end-to-end --very-sensitive --no-mixed --no-discordant --phred33 -I 10 -X 700 -p ${cores} -x ${spikeInRef} -1 ${projPath}/fastq/${histName}_R1.fastq.gz -2 ${projPath}/fastq/${histName}_R2.fastq.gz -S $projPath/alignment/sam/${histName}_bowtie2_spikeIn.sam &amp;&gt; $projPath/alignment/sam/bowtie2_summary/${histName}_bowtie2_spikeIn.txt

seqDepthDouble=`samtools view -F 0x04 $projPath/alignment/sam/${histName}_bowtie2_spikeIn.sam | wc -l`
seqDepth=$((seqDepthDouble/2))
echo $seqDepth &gt;$projPath/alignment/sam/bowtie2_summary/${histName}_bowtie2_spikeIn.seqDepth
</code></pre>
<ul>
<li>For spike-in normalization, reads are aligned to the E. coli genome U00096.3 with two more parameters <code>--no-overlap</code> and <code>--no-dovetail</code> (<code>--end-to-end --very-sensitive --no-overlap --no-dovetail --no-mixed --no-discordant --phred33 -I 10 -X 700</code>) to avoid possible cross-mapping of the experimental genome to that of the carry-over E. coli DNA that is used for calibration.</li>
</ul>
</div>
<div id="alignment-summary" class="section level3">
<h3>3.1.3 Alignment summary</h3>
<p>For more detailed parameters explanation, users can refer to the <a href="http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml">bowite2 manual</a>.</p>
<p>Bowtie2 alignment results summary is saved at <code>${projPath}/alignment/sam/bowtie2)summary/${histName}_bowtie2.txt</code> and you should expect the results look similar.</p>
<pre><code>2984630 reads; of these:
  2984630 (100.00%) were paired; of these:
    125110 (4.19%) aligned concordantly 0 times
    2360430 (79.09%) aligned concordantly exactly 1 time
    499090 (16.72%) aligned concordantly &gt;1 times
95.81% overall alignment rate</code></pre>
<ul>
<li>2984640 is the sequencing depth, i.e., total number of paired reads.</li>
<li>125110 is the number of read pairs that fail to be mapped.</li>
<li>2360430 + 499090 is the number of read paris that are successfully mapped.</li>
<li>95.81% is the overall alignment rate</li>
</ul>
</div>
</div>
<div id="report-sequencing-mapping-summary-required" class="section level2">
<h2>3.2 Report sequencing mapping summary [required]</h2>
<p>Summarize the raw reads and uniquely mapping reads to report the efficiency of alignment. Alignment frequencies are expected to be &gt;80% for high-quality data. CUT&amp;Tag data typically has very low backgrounds, so as few as 1 million mapped fragments can give robust profiles for a histone modification in the human genome. Profiling of less-abundant transcription factors and chromatin proteins may require 10 times as many mapped fragments for downstream analysis.</p>
<p>We can evaluate the following metrics:</p>
<ul>
<li>Sequencing depth</li>
<li>Alignment rate</li>
<li>Number of mappable fragments</li>
<li>Duplication rate</li>
<li>Unique library size</li>
<li>Fragment size distribution</li>
</ul>
<div id="sequencing-depth" class="section level3">
<h3>3.2.1 Sequencing depth</h3>
<pre class="r"><code>##=== R command ===## 
## Path to the project and histone list
projPath = &quot;/fh/fast/gottardo_r/yezheng_working/cuttag/CUTTag_tutorial&quot;
sampleList = c(&quot;K27me3_rep1&quot;, &quot;K27me3_rep2&quot;, &quot;K4me3_rep1&quot;, &quot;K4me3_rep2&quot;, &quot;IgG_rep1&quot;, &quot;IgG_rep2&quot;)
histList = c(&quot;K27me3&quot;, &quot;K4me3&quot;, &quot;IgG&quot;)

## Collect the alignment results from the bowtie2 alignment summary files
alignResult = c()
for(hist in sampleList){
  alignRes = read.table(paste0(projPath, &quot;/alignment/sam/bowtie2_summary/&quot;, hist, &quot;_bowtie2.txt&quot;), header = FALSE, fill = TRUE)
  alignRate = substr(alignRes$V1[6], 1, nchar(as.character(alignRes$V1[6]))-1)
  histInfo = strsplit(hist, &quot;_&quot;)[[1]]
  alignResult = data.frame(Histone = histInfo[1], Replicate = histInfo[2], 
                           SequencingDepth = alignRes$V1[1] %&gt;% as.character %&gt;% as.numeric, 
                           MappedFragNum_hg38 = alignRes$V1[4] %&gt;% as.character %&gt;% as.numeric + alignRes$V1[5] %&gt;% as.character %&gt;% as.numeric, 
                           AlignmentRate_hg38 = alignRate %&gt;% as.numeric)  %&gt;% rbind(alignResult, .)
}
alignResult$Histone = factor(alignResult$Histone, levels = histList)
alignResult %&gt;% mutate(AlignmentRate_hg38 = paste0(AlignmentRate_hg38, &quot;%&quot;))</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["fctr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["SequencingDepth"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["MappedFragNum_hg38"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_hg38"],"name":[5],"type":["chr"],"align":["left"]}],"data":[{"1":"K27me3","2":"rep1","3":"2984630","4":"2859520","5":"95.81%"},{"1":"K27me3","2":"rep2","3":"2702260","4":"2606295","5":"96.45%"},{"1":"K4me3","2":"rep1","3":"1581710","4":"1494122","5":"94.46%"},{"1":"K4me3","2":"rep2","3":"1885056","4":"1742005","5":"92.41%"},{"1":"IgG","2":"rep1","3":"2127635","4":"1747065","5":"82.11%"},{"1":"IgG","2":"rep2","3":"2192908","4":"1992929","5":"90.88%"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
</div>
<div id="spike-in-alignment" class="section level3">
<h3>3.2.2 Spike-in alignment</h3>
<pre class="r"><code>##=== R command ===## 
spikeAlign = c()
for(hist in sampleList){
  spikeRes = read.table(paste0(projPath, &quot;/alignment/sam/bowtie2_summary/&quot;, hist, &quot;_bowtie2_spikeIn.txt&quot;), header = FALSE, fill = TRUE)
  alignRate = substr(spikeRes$V1[6], 1, nchar(as.character(spikeRes$V1[6]))-1)
  histInfo = strsplit(hist, &quot;_&quot;)[[1]]
  spikeAlign = data.frame(Histone = histInfo[1], Replicate = histInfo[2], 
                          SequencingDepth = spikeRes$V1[1] %&gt;% as.character %&gt;% as.numeric, 
                          MappedFragNum_spikeIn = spikeRes$V1[4] %&gt;% as.character %&gt;% as.numeric + spikeRes$V1[5] %&gt;% as.character %&gt;% as.numeric, 
                          AlignmentRate_spikeIn = alignRate %&gt;% as.numeric)  %&gt;% rbind(spikeAlign, .)
}
spikeAlign$Histone = factor(spikeAlign$Histone, levels = histList)
spikeAlign %&gt;% mutate(AlignmentRate_spikeIn = paste0(AlignmentRate_spikeIn, &quot;%&quot;))</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["fctr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["SequencingDepth"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["MappedFragNum_spikeIn"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_spikeIn"],"name":[5],"type":["chr"],"align":["left"]}],"data":[{"1":"K27me3","2":"rep1","3":"2984630","4":"235","5":"0.01%"},{"1":"K27me3","2":"rep2","3":"2702260","4":"487","5":"0.02%"},{"1":"K4me3","2":"rep1","3":"1581710","4":"375","5":"0.02%"},{"1":"K4me3","2":"rep2","3":"1885056","4":"4442","5":"0.24%"},{"1":"IgG","2":"rep1","3":"2127635","4":"75733","5":"3.56%"},{"1":"IgG","2":"rep2","3":"2192908","4":"79123","5":"3.61%"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
</div>
<div id="summarize-the-alignment-to-hg38-and-e.coli" class="section level3">
<h3>3.2.3 Summarize the alignment to hg38 and E.coli</h3>
<pre class="r"><code>##=== R command ===## 
alignSummary = left_join(alignResult, spikeAlign, by = c(&quot;Histone&quot;, &quot;Replicate&quot;, &quot;SequencingDepth&quot;)) %&gt;%
  mutate(AlignmentRate_hg38 = paste0(AlignmentRate_hg38, &quot;%&quot;), 
         AlignmentRate_spikeIn = paste0(AlignmentRate_spikeIn, &quot;%&quot;))
alignSummary</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["fctr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["SequencingDepth"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["MappedFragNum_hg38"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_hg38"],"name":[5],"type":["chr"],"align":["left"]},{"label":["MappedFragNum_spikeIn"],"name":[6],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_spikeIn"],"name":[7],"type":["chr"],"align":["left"]}],"data":[{"1":"K27me3","2":"rep1","3":"2984630","4":"2859520","5":"95.81%","6":"235","7":"0.01%"},{"1":"K27me3","2":"rep2","3":"2702260","4":"2606295","5":"96.45%","6":"487","7":"0.02%"},{"1":"K4me3","2":"rep1","3":"1581710","4":"1494122","5":"94.46%","6":"375","7":"0.02%"},{"1":"K4me3","2":"rep2","3":"1885056","4":"1742005","5":"92.41%","6":"4442","7":"0.24%"},{"1":"IgG","2":"rep1","3":"2127635","4":"1747065","5":"82.11%","6":"75733","7":"3.56%"},{"1":"IgG","2":"rep2","3":"2192908","4":"1992929","5":"90.88%","6":"79123","7":"3.61%"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
</div>
<div id="visualizing-the-sequencing-depth-and-alignment-results." class="section level3">
<h3>3.2.4 Visualizing the sequencing depth and alignment results.</h3>
<pre class="r"><code>##=== R command ===## 
## Generate sequencing depth boxplot
fig3A = alignResult %&gt;% ggplot(aes(x = Histone, y = SequencingDepth/1000000, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;Sequencing Depth per Million&quot;) +
    xlab(&quot;&quot;) + 
    ggtitle(&quot;A. Sequencing Depth&quot;)

fig3B = alignResult %&gt;% ggplot(aes(x = Histone, y = MappedFragNum_hg38/1000000, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;Mapped Fragments per Million&quot;) +
    xlab(&quot;&quot;) +
    ggtitle(&quot;B. Alignable Fragment (hg38)&quot;)

fig3C = alignResult %&gt;% ggplot(aes(x = Histone, y = AlignmentRate_hg38, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;% of Mapped Fragments&quot;) +
    xlab(&quot;&quot;) +
    ggtitle(&quot;C. Alignment Rate (hg38)&quot;)

fig3D = spikeAlign %&gt;% ggplot(aes(x = Histone, y = AlignmentRate_spikeIn, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;Spike-in Alignment Rate&quot;) +
    xlab(&quot;&quot;) +
    ggtitle(&quot;D. Alignment Rate (E.coli)&quot;)

ggarrange(fig3A, fig3B, fig3C, fig3D, ncol = 2, nrow=2, common.legend = TRUE, legend=&quot;bottom&quot;)</code></pre>
<p><img src="figure/tutorials.Rmd/unnamed-chunk-12-1.png" width="100%" style="display: block; margin: auto;" /></p>
<p>
<button type="button" class="btn btn-default btn-xs btn-workflowr btn-workflowr-fig" data-toggle="collapse" data-target="#fig-unnamed-chunk-12-1">
Past versions of unnamed-chunk-12-1.png
</button>
</p>
<div id="fig-unnamed-chunk-12-1" class="collapse">
<div class="table-responsive">
<table class="table table-condensed table-hover">
<thead>
<tr>
<th>
Version
</th>
<th>
Author
</th>
<th>
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/b2993411ac609462d4ebc4b7ce1446e327a00e56/docs/figure/tutorials.Rmd/unnamed-chunk-12-1.png" target="_blank">b299341</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/d7f438eb1c959c4b627814b1aa1f07850756b6bd/docs/figure/tutorials.Rmd/unnamed-chunk-12-1.png" target="_blank">d7f438e</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
</tbody>
</table>
</div>
</div>
<p>In a typical CUT&amp;Tag experiment targeting the abundant H3K27me3 histone modification in 65,000 K562 cells, the percentage of E. coli reads range from ~0.01% to 10%. With fewer cells or less abundant epitopes, E. coli reads can comprise as much as 70% or the total mapped reads. For IgG controls, the percentage of E. coli reads is typically much higher than that for an abundant histone modification.</p>
</div>
</div>
<div id="remove-duplicates-optionalrequired" class="section level2">
<h2>3.3. Remove duplicates [optional/required]</h2>
<p>CUT&amp;Tag integrates adapters into DNA in the vicinity of the antibody-tethered pA-Tn5, and the exact sites of integration are affected by the accessibility of surrounding DNA. For this reason fragments that share exact starting and ending positions are expected to be common, and such ‘duplicates’ may not be due to duplication during PCR. In practice, we have found that the apparent duplication rate is low for high quality CUT&amp;Tag datasets, and even the apparent ‘duplicate’ fragments are likely to be true fragments. Thus, we do not recommend removing the duplicates. In experiments with very small amounts of material or where PCR duplication is suspected, duplicates can be removed. The following commands show how to check the duplication rate using <a href="https://broadinstitute.github.io/picard/">Picard</a>.</p>
<pre class="bash"><code>##== linux command ==##
## depending on how you load picard and your server environment, the picardCMD can be different. Adjust accordingly.
picardCMD=&quot;java -jar picard.jar&quot;
mkdir -p $projPath/alignment/removeDuplicate/picard_summary

## Sort by coordinate
$picardCMD SortSam I=$projPath/alignment/sam/${histName}_bowtie2.sam O=$projPath/alignment/sam/${histName}_bowtie2.sorted.sam SORT_ORDER=coordinate

## mark duplicates
$picardCMD MarkDuplicates I=$projPath/alignment/sam/${histName}_bowtie2.sorted.sam O=$projPath/alignment/removeDuplicate/${histName}_bowtie2.sorted.dupMarked.sam METRICS_FILE=$projPath/alignment/removeDuplicate/picard_summary/${histName}_picard.dupMark.txt

## remove duplicates
picardCMD MarkDuplicates I=$projPath/alignment/sam/${histName}_bowtie2.sorted.sam O=$projPath/alignment/removeDuplicate/${histName}_bowtie2.sorted.rmDup.sam REMOVE_DUPLICATES=true METRICS_FILE=$projPath/alignment/removeDuplicate/picard_summary/${histName}_picard.rmDup.txt</code></pre>
<p>We summarize the apparent duplication rate and calculate the unique library size without duplicates.</p>
<pre class="r"><code>##=== R command ===## 
## Summarize the duplication information from the picard summary outputs.
dupResult = c()
for(hist in sampleList){
  dupRes = read.table(paste0(projPath, &quot;/alignment/removeDuplicate/picard_summary/&quot;, hist, &quot;_picard.rmDup.txt&quot;), header = TRUE, fill = TRUE)
  
  histInfo = strsplit(hist, &quot;_&quot;)[[1]]
  dupResult = data.frame(Histone = histInfo[1], Replicate = histInfo[2], MappedFragNum_hg38 = dupRes$READ_PAIRS_EXAMINED[1] %&gt;% as.character %&gt;% as.numeric, DuplicationRate = dupRes$PERCENT_DUPLICATION[1] %&gt;% as.character %&gt;% as.numeric * 100, EstimatedLibrarySize = dupRes$ESTIMATED_LIBRARY_SIZE[1] %&gt;% as.character %&gt;% as.numeric) %&gt;% mutate(UniqueFragNum = MappedFragNum_hg38 * (1-DuplicationRate/100))  %&gt;% rbind(dupResult, .)
}
dupResult$Histone = factor(dupResult$Histone, levels = histList)
alignDupSummary = left_join(alignSummary, dupResult, by = c(&quot;Histone&quot;, &quot;Replicate&quot;, &quot;MappedFragNum_hg38&quot;)) %&gt;% mutate(DuplicationRate = paste0(DuplicationRate, &quot;%&quot;))
alignDupSummary</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["fctr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["SequencingDepth"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["MappedFragNum_hg38"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_hg38"],"name":[5],"type":["chr"],"align":["left"]},{"label":["MappedFragNum_spikeIn"],"name":[6],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_spikeIn"],"name":[7],"type":["chr"],"align":["left"]},{"label":["DuplicationRate"],"name":[8],"type":["chr"],"align":["left"]},{"label":["EstimatedLibrarySize"],"name":[9],"type":["dbl"],"align":["right"]},{"label":["UniqueFragNum"],"name":[10],"type":["dbl"],"align":["right"]}],"data":[{"1":"K27me3","2":"rep1","3":"2984630","4":"2859520","5":"95.81%","6":"235","7":"0.01%","8":"4.8503%","9":"28683768","10":"2720824.7"},{"1":"K27me3","2":"rep2","3":"2702260","4":"2606295","5":"96.45%","6":"487","7":"0.02%","8":"1.0445%","9":"125234408","10":"2579072.2"},{"1":"K4me3","2":"rep1","3":"1581710","4":"1494122","5":"94.46%","6":"375","7":"0.02%","8":"6.6056%","9":"10843753","10":"1395426.3"},{"1":"K4me3","2":"rep2","3":"1885056","4":"1742005","5":"92.41%","6":"4442","7":"0.24%","8":"2.7026%","9":"31790431","10":"1694925.6"},{"1":"IgG","2":"rep1","3":"2127635","4":"1747065","5":"82.11%","6":"75733","7":"3.56%","8":"81.336%","9":"327661","10":"326072.2"},{"1":"IgG","2":"rep2","3":"2192908","4":"1992929","5":"90.88%","6":"79123","7":"3.61%","8":"34.3392%","9":"2192721","10":"1308573.1"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
<ul>
<li><p>In these example datasets, the IgG control samples have relatively high duplication rates, since reads in this sample derive from non-specific tagmentation in the CUT&amp;Tag reactions. Therefore, it is appropriate to remove the duplicates from the IgG datasets before downstream analysis.</p></li>
<li><p>The estimated library size are the estimated number of unique molecules in the library based on PE duplication calculated by Picard.</p></li>
<li><p>The estimated library sizes is proportional to the abundance of the targeted epitope and to the quality of the antibody used, while the estimated library sizes of IgG samples are expected to be very low.</p></li>
<li><p>Unique fragment number is calculated by the MappedFragNum_hg38 * (1-DuplicationRate/100).</p></li>
</ul>
<pre class="r"><code>##=== R command ===## 
## generate boxplot figure for the  duplication rate
fig4A = dupResult %&gt;% ggplot(aes(x = Histone, y = DuplicationRate, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;Duplication Rate (*100%)&quot;) +
    xlab(&quot;&quot;) 

fig4B = dupResult %&gt;% ggplot(aes(x = Histone, y = EstimatedLibrarySize, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;Estimated Library Size&quot;) +
    xlab(&quot;&quot;) 

fig4C = dupResult %&gt;% ggplot(aes(x = Histone, y = UniqueFragNum, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;# of Unique Fragments&quot;) +
    xlab(&quot;&quot;)

ggarrange(fig4A, fig4B, fig4C, ncol = 3, common.legend = TRUE, legend=&quot;bottom&quot;)</code></pre>
<p><img src="figure/tutorials.Rmd/unnamed-chunk-15-1.png" width="100%" style="display: block; margin: auto;" /></p>
<p>
<button type="button" class="btn btn-default btn-xs btn-workflowr btn-workflowr-fig" data-toggle="collapse" data-target="#fig-unnamed-chunk-15-1">
Past versions of unnamed-chunk-15-1.png
</button>
</p>
<div id="fig-unnamed-chunk-15-1" class="collapse">
<div class="table-responsive">
<table class="table table-condensed table-hover">
<thead>
<tr>
<th>
Version
</th>
<th>
Author
</th>
<th>
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/b2993411ac609462d4ebc4b7ce1446e327a00e56/docs/figure/tutorials.Rmd/unnamed-chunk-15-1.png" target="_blank">b299341</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/d7f438eb1c959c4b627814b1aa1f07850756b6bd/docs/figure/tutorials.Rmd/unnamed-chunk-15-1.png" target="_blank">d7f438e</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
</tbody>
</table>
</div>
</div>
</div>
<div id="assess-mapped-fragment-size-distribution-required" class="section level2">
<h2>3.4. Assess mapped fragment size distribution [Required]</h2>
<p>CUT&amp;Tag inserts adapters on either side of chromatin particles in the vicinity of the tethered enzyme, although tagmentation within chromatin particles can also occur. So, CUT&amp;Tag reactions targeting a histone modification predominantly results in fragments that are nucleosomal lengths (~180 bp), or multiples of that length. CUT&amp;Tag targeting transcription factors predominantly produce nucleosome-sized fragments and variable amounts of shorter fragments, from neighboring nucleosomes and the factor-bound site, respectively. Tagmentation of DNA on the surface of nucleosomes also occurs, and plotting fragment lengths with single-basepair resolution reveal a 10-bp sawtooth periodicity, which is typical of successful CUT&amp;Tag experiments.</p>
<pre class="bash"><code>##== linux command ==##
mkdir -p $projPath/alignment/sam/fragmentLen

## Extract the 9th column from the alignment sam file which is the fragment length
samtools view -F 0x04 $projPath/alignment/sam/${histName}_bowtie2.sam | awk -F&#39;\t&#39; &#39;function abs(x){return ((x &lt; 0.0) ? -x : x)} {print abs($9)}&#39; | sort | uniq -c | awk -v OFS=&quot;\t&quot; &#39;{print $2, $1/2}&#39; &gt;$projPath/alignment/sam/fragmentLen/${histName}_fragmentLen.txt</code></pre>
<pre class="r"><code>##=== R command ===## 
## Collect the fragment size information
fragLen = c()
for(hist in sampleList){
  
  histInfo = strsplit(hist, &quot;_&quot;)[[1]]
  fragLen = read.table(paste0(projPath, &quot;/alignment/sam/fragmentLen/&quot;, hist, &quot;_fragmentLen.txt&quot;), header = FALSE) %&gt;% mutate(fragLen = V1 %&gt;% as.numeric, fragCount = V2 %&gt;% as.numeric, Weight = as.numeric(V2)/sum(as.numeric(V2)), Histone = histInfo[1], Replicate = histInfo[2], sampleInfo = hist) %&gt;% rbind(fragLen, .) 
}
fragLen$sampleInfo = factor(fragLen$sampleInfo, levels = sampleList)
fragLen$Histone = factor(fragLen$Histone, levels = histList)
## Generate the fragment size density plot (violin plot)
fig5A = fragLen %&gt;% ggplot(aes(x = sampleInfo, y = fragLen, weight = Weight, fill = Histone)) +
    geom_violin(bw = 5) +
    scale_y_continuous(breaks = seq(0, 800, 50)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 20) +
    ggpubr::rotate_x_text(angle = 20) +
    ylab(&quot;Fragment Length&quot;) +
    xlab(&quot;&quot;)

fig5B = fragLen %&gt;% ggplot(aes(x = fragLen, y = fragCount, color = Histone, group = sampleInfo, linetype = Replicate)) +
  geom_line(size = 1) +
  scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;) +
  theme_bw(base_size = 20) +
  xlab(&quot;Fragment Length&quot;) +
  ylab(&quot;Count&quot;) +
  coord_cartesian(xlim = c(0, 500))

ggarrange(fig5A, fig5B, ncol = 2)</code></pre>
<p><img src="figure/tutorials.Rmd/unnamed-chunk-17-1.png" width="100%" style="display: block; margin: auto;" /></p>
<p>
<button type="button" class="btn btn-default btn-xs btn-workflowr btn-workflowr-fig" data-toggle="collapse" data-target="#fig-unnamed-chunk-17-1">
Past versions of unnamed-chunk-17-1.png
</button>
</p>
<div id="fig-unnamed-chunk-17-1" class="collapse">
<div class="table-responsive">
<table class="table table-condensed table-hover">
<thead>
<tr>
<th>
Version
</th>
<th>
Author
</th>
<th>
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/b2993411ac609462d4ebc4b7ce1446e327a00e56/docs/figure/tutorials.Rmd/unnamed-chunk-17-1.png" target="_blank">b299341</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/d7f438eb1c959c4b627814b1aa1f07850756b6bd/docs/figure/tutorials.Rmd/unnamed-chunk-17-1.png" target="_blank">d7f438e</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
</tbody>
</table>
</div>
</div>
<ul>
<li>The smaller fragments (50-100 bp) can be due to that tethered Tn5 can tagment on the surface of a nucleosome as well as in linker regions, so the small fragments might not be background.</li>
</ul>
</div>
<div id="assess-replicate-reproducibility" class="section level2">
<h2>3.5 Assess replicate reproducibility</h2>
<p>Data reproducibility between replicates is assessed by correlation analysis of mapped read counts across the genome. For the simplicity of implementation, we will postpone this analysis after Section IV when the file format has been converted into fragment bed files.</p>
</div>
</div>
<div id="iv.-alignment-results-filtering-and-file-format-conversion" class="section level1">
<h1>IV. Alignment results filtering and file format conversion</h1>
<div id="filtering-mapped-reads-by-the-mapping-quality-filtering-optinal" class="section level2">
<h2>4.1 Filtering mapped reads by the mapping quality filtering [optinal]</h2>
<p>Some project may require more stringent filtering on the alignment quality score. This <a href="http://biofinysics.blogspot.com/2014/05/how-does-bowtie2-assign-mapq-scores.html">blog</a> detailedly discussed how does bowtie assign quality score with examples.</p>
<p>MAPQ(x) = -10 * <span class="math inline">\(log_{10}\)</span>(P(x is mapped wrongly)) = -10 * <span class="math inline">\(log_{10}(p)\)</span></p>
<p>which ranges from 0 to 37, 40 or 42.</p>
<p><code>samtools view -q minQualityScore</code> will eliminate all the alignment results that are below the minQualityScore defined by user.</p>
<pre class="bash"><code>##== linux command ==##
minQualityScore=2
samtools view -q $minQualityScore ${projPath}/alignment/sam/${histName}_bowtie2.sam &gt;${projPath}/alignment/sam/${histName}_bowtie2.qualityScore$minQualityScore.sam</code></pre>
<ul>
<li>If you do implement this filtering, please replace the <code>${histName}_bowtie2.sam</code> in the following steps by this filtered sam file <code>${histName}_bowtie2.qualityScore$minQualityScore.sam</code>.</li>
</ul>
</div>
<div id="file-format-conversion-required" class="section level2">
<h2>4.2 File format conversion [required]</h2>
<p>This section is <strong>required</strong> in preparation for the peak calling and visualization where there are a few filtering and file format conversion that need to be done.</p>
<pre class="bash"><code>##== linux command ==##
## Filter and keep the mapped read pairs
samtools view -bS -F 0x04 $projPath/alignment/sam/${histName}_bowtie2.sam &gt;$projPath/alignment/bam/${histName}_bowtie2.mapped.bam

## Convert into bed file format
bedtools bamtobed -i $projPath/alignment/bam/${histName}_bowtie2.mapped.bam -bedpe &gt;$projPath/alignment/bed/${histName}_bowtie2.bed

## Keep the read pairs that are on the same chromosome and fragment length less than 1000bp.
awk &#39;$1==$4 &amp;&amp; $6-$2 &lt; 1000 {print $0}&#39; $projPath/alignment/bed/${histName}_bowtie2.bed &gt;$projPath/alignment/bed/${histName}_bowtie2.clean.bed

## Only extract the fragment related columns
cut -f 1,2,6 $projPath/alignment/bed/${histName}_bowtie2.clean.bed | sort -k1,1 -k2,2n -k3,3n  &gt;$projPath/alignment/bed/${histName}_bowtie2.fragments.bed</code></pre>
</div>
<div id="assess-replicate-reproducibility-continue-section-3.5" class="section level2">
<h2>4.3 Assess replicate reproducibility (continue section 3.5)</h2>
<p>To study the reproducibility between replicates and across conditions, the genome is split into 500 bp bins, and a Pearson correlation of the log2-transformed values of read counts in each bin is calculated between replicate datasets. Multiple replicates and IgG control datasets are displayed in a hierarchically clustered correlation matrix.</p>
<pre class="bash"><code>##== linux command ==##
## We use the mid point of each fragment to infer which 500bp bins does this fragment belong to.
binLen=500
awk -v w=$binLen &#39;{print $1, int(($2 + $3)/(2*w))*w + w/2}&#39; $projPath/alignment/bed/${histName}_bowtie2.fragments.bed | sort -k1,1V -k2,2n | uniq -c | awk -v OFS=&quot;\t&quot; &#39;{print $2, $3, $1}&#39; |  sort -k1,1V -k2,2n  &gt;$projPath/alignment/bed/${histName}_bowtie2.fragmentsCount.bin$binLen.bed
</code></pre>
<pre class="r"><code>##== R command ==##
reprod = c()
fragCount = NULL
for(hist in sampleList){
  
  if(is.null(fragCount)){
    
    fragCount = read.table(paste0(projPath, &quot;/alignment/bed/&quot;, hist, &quot;_bowtie2.fragmentsCount.bin500.bed&quot;), header = FALSE) 
    colnames(fragCount) = c(&quot;chrom&quot;, &quot;bin&quot;, hist)
  
  }else{
    
    fragCountTmp = read.table(paste0(projPath, &quot;/alignment/bed/&quot;, hist, &quot;_bowtie2.fragmentsCount.bin500.bed&quot;), header = FALSE)
    colnames(fragCountTmp) = c(&quot;chrom&quot;, &quot;bin&quot;, hist)
    fragCount = full_join(fragCount, fragCountTmp, by = c(&quot;chrom&quot;, &quot;bin&quot;))
    
  }
}

M = cor(fragCount %&gt;% select(-c(&quot;chrom&quot;, &quot;bin&quot;)) %&gt;% log2(), use = &quot;complete.obs&quot;) 

corrplot(M, method = &quot;color&quot;, outline = T, addgrid.col = &quot;darkgray&quot;, order=&quot;hclust&quot;, addrect = 3, rect.col = &quot;black&quot;, rect.lwd = 3,cl.pos = &quot;b&quot;, tl.col = &quot;indianred4&quot;, tl.cex = 1, cl.cex = 1, addCoef.col = &quot;black&quot;, number.digits = 2, number.cex = 1, col = colorRampPalette(c(&quot;midnightblue&quot;,&quot;white&quot;,&quot;darkred&quot;))(100))</code></pre>
<p><img src="figure/tutorials.Rmd/unnamed-chunk-21-1.png" width="100%" style="display: block; margin: auto;" /></p>
<p>
<button type="button" class="btn btn-default btn-xs btn-workflowr btn-workflowr-fig" data-toggle="collapse" data-target="#fig-unnamed-chunk-21-1">
Past versions of unnamed-chunk-21-1.png
</button>
</p>
<div id="fig-unnamed-chunk-21-1" class="collapse">
<div class="table-responsive">
<table class="table table-condensed table-hover">
<thead>
<tr>
<th>
Version
</th>
<th>
Author
</th>
<th>
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/b2993411ac609462d4ebc4b7ce1446e327a00e56/docs/figure/tutorials.Rmd/unnamed-chunk-21-1.png" target="_blank">b299341</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/d544756d9ed67d5e7ab554eb1fd2fc9e410c5f13/docs/figure/tutorials.Rmd/unnamed-chunk-21-1.png" target="_blank">d544756</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-06-01
</td>
</tr>
</tbody>
</table>
</div>
</div>
</div>
</div>
<div id="v.-spike-in-calibration" class="section level1">
<h1>V. Spike-in calibration</h1>
<p>This section is <strong>optional</strong> but <strong>recommended</strong> depending on your experimental protocol. We have shown the alignment to the spike-in genome in Section 3.1.2 and the spike-in alignment summary in Section 3.2.2.</p>
<p>The underlying assumption is that the ratio of fragments mapped to the primary genome to the E. coli genome is the same for a series of samples, each using the same number of cells. Because of this assumption, we do not normalize between experiments or between batches of purified pATn5, which can have very different amounts of carry-over E. coli DNA. Using a constant C to avoid small fractions in normalized data, we define a scaling factor S as</p>
<p><code>S = C / (fragments mapped to E. coli genome)</code></p>
<p>Normalized coverage is then calculated as:</p>
<p><code>Normalized coverage = (primary_genome_coverage) * S</code></p>
<p>The Constant is an arbitrary multiplier, typically 10,000. The resulting file will be comparatively small as a genomic coverage bedgraph file.</p>
<pre class="bash"><code>##== linux command ==##
if [[ &quot;$seqDepth&quot; -gt &quot;1&quot; ]]; then
    
    mkdir -p $projPath/alignment/bedgraph

    scale_factor=`echo &quot;10000 / $seqDepth&quot; | bc -l`
    echo &quot;Scaling factor for $histName is: $scale_factor!&quot;
    bedtools genomecov -bg -scale $scale_factor -i $projPath/alignment/bed/${histName}_bowtie2.fragments.bed -g $chromSize &gt; $projPath/alignment/bedgraph/${histName}_bowtie2.fragments.normalized.bedgraph
    
fi</code></pre>
<div id="scaling-factor" class="section level2">
<h2>5.1 Scaling factor</h2>
<pre class="r"><code>##=== R command ===## 
scaleFactor = c()
multiplier = 10000
for(hist in sampleList){
  spikeDepth = read.table(paste0(projPath, &quot;/alignment/sam/bowtie2_summary/&quot;, hist, &quot;_bowtie2_spikeIn.seqDepth&quot;), header = FALSE, fill = TRUE)$V1[1]
  
  histInfo = strsplit(hist, &quot;_&quot;)[[1]]
  scaleFactor = data.frame(scaleFactor = multiplier/spikeDepth, Histone = histInfo[1], Replicate = histInfo[2])  %&gt;% rbind(scaleFactor, .)
}
scaleFactor$Histone = factor(scaleFactor$Histone, levels = histList)
left_join(alignDupSummary, scaleFactor, by = c(&quot;Histone&quot;, &quot;Replicate&quot;))</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["fctr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["SequencingDepth"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["MappedFragNum_hg38"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_hg38"],"name":[5],"type":["chr"],"align":["left"]},{"label":["MappedFragNum_spikeIn"],"name":[6],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_spikeIn"],"name":[7],"type":["chr"],"align":["left"]},{"label":["DuplicationRate"],"name":[8],"type":["chr"],"align":["left"]},{"label":["EstimatedLibrarySize"],"name":[9],"type":["dbl"],"align":["right"]},{"label":["UniqueFragNum"],"name":[10],"type":["dbl"],"align":["right"]},{"label":["scaleFactor"],"name":[11],"type":["dbl"],"align":["right"]}],"data":[{"1":"K27me3","2":"rep1","3":"2984630","4":"2859520","5":"95.81%","6":"235","7":"0.01%","8":"4.8503%","9":"28683768","10":"2720824.7","11":"42.5531915"},{"1":"K27me3","2":"rep2","3":"2702260","4":"2606295","5":"96.45%","6":"487","7":"0.02%","8":"1.0445%","9":"125234408","10":"2579072.2","11":"20.5338809"},{"1":"K4me3","2":"rep1","3":"1581710","4":"1494122","5":"94.46%","6":"375","7":"0.02%","8":"6.6056%","9":"10843753","10":"1395426.3","11":"26.6666667"},{"1":"K4me3","2":"rep2","3":"1885056","4":"1742005","5":"92.41%","6":"4442","7":"0.24%","8":"2.7026%","9":"31790431","10":"1694925.6","11":"2.2512382"},{"1":"IgG","2":"rep1","3":"2127635","4":"1747065","5":"82.11%","6":"75733","7":"3.56%","8":"81.336%","9":"327661","10":"326072.2","11":"0.1320428"},{"1":"IgG","2":"rep2","3":"2192908","4":"1992929","5":"90.88%","6":"79123","7":"3.61%","8":"34.3392%","9":"2192721","10":"1308573.1","11":"0.1263855"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
<pre class="r"><code>##=== R command ===##
## Generate sequencing depth boxplot
fig6A = scaleFactor %&gt;% ggplot(aes(x = Histone, y = scaleFactor, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 20) +
    ylab(&quot;Spike-in Scalling Factor&quot;) +
    xlab(&quot;&quot;)

normDepth = inner_join(scaleFactor, alignResult, by = c(&quot;Histone&quot;, &quot;Replicate&quot;)) %&gt;% mutate(normDepth = MappedFragNum_hg38 * scaleFactor)

fig6B = normDepth %&gt;% ggplot(aes(x = Histone, y = normDepth, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.9, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 20) +
    ylab(&quot;Normalization Fragment Count&quot;) +
    xlab(&quot;&quot;) + 
    coord_cartesian(ylim = c(1000000, 130000000))
ggarrange(fig6A, fig6B, ncol = 2, common.legend = TRUE, legend=&quot;bottom&quot;)</code></pre>
<p><img src="figure/tutorials.Rmd/unnamed-chunk-24-1.png" width="100%" style="display: block; margin: auto;" /></p>
<p>
<button type="button" class="btn btn-default btn-xs btn-workflowr btn-workflowr-fig" data-toggle="collapse" data-target="#fig-unnamed-chunk-24-1">
Past versions of unnamed-chunk-24-1.png
</button>
</p>
<div id="fig-unnamed-chunk-24-1" class="collapse">
<div class="table-responsive">
<table class="table table-condensed table-hover">
<thead>
<tr>
<th>
Version
</th>
<th>
Author
</th>
<th>
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/b2993411ac609462d4ebc4b7ce1446e327a00e56/docs/figure/tutorials.Rmd/unnamed-chunk-24-1.png" target="_blank">b299341</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/d544756d9ed67d5e7ab554eb1fd2fc9e410c5f13/docs/figure/tutorials.Rmd/unnamed-chunk-24-1.png" target="_blank">d544756</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-06-01
</td>
</tr>
</tbody>
</table>
</div>
</div>
</div>
</div>
<div id="vi.-peak-calling" class="section level1">
<h1>VI. Peak calling</h1>
<div id="seacr" class="section level2">
<h2>6.1. SEACR</h2>
<p>The Sparse Enrichment Analysis for CUT&amp;RUN, <a href="https://github.com/FredHutch/SEACR/">SEACR</a>, package is designed to call peaks and enriched regions from chromatin profiling data with very low backgrounds (i.e., regions with no read coverage) that are typical for CUT&amp;Tag chromatin profiling experiments. SEACR requires bedGraph files from paired-end sequencing as input and defines peaks as contiguous blocks of basepair coverage that do not overlap with blocks of background signal delineated in the IgG control dataset. SEACR is effective for calling both narrow peaks from factor binding sites and broad domains characteristic of some histone modifications. The description of the method is published at <a href="https://epigeneticsandchromatin.biomedcentral.com/articles/10.1186/s13072-019-0287-4">Meers et al. 2019</a>, and the user’s manual is available on <a href="https://github.com/FredHutch/SEACR/">github</a>. Since we have normalized fragment counts with the E. coli read count, we set the normalization option of SEACR to “non”. Otherwise, the “norm” is recommended.</p>
<pre class="bash"><code>##== linux command ==##
seacr=&quot;/fh/fast/gottardo_r/yezheng_working/Software/SEACR/SEACR_1.3.sh&quot;
histControl=$2
mkdir -p $projPath/peakCalling/SEACR

bash $seacr $projPath/alignment/bedgraph/${histName}_bowtie2.fragments.normalized.bedgraph \
     $projPath/alignment/bedgraph/${histControl}_bowtie2.fragments.normalized.bedgraph \
     non stringent $projPath/peakCalling/SEACR/${histName}_seacr_control.peaks

bash $seacr $projPath/alignment/bedgraph/${histName}_bowtie2.fragments.normalized.bedgraph 0.01 non stringent $projPath/peakCalling/SEACR/${histName}_seacr_top0.01.peaks</code></pre>
<div id="number-of-peaks-called" class="section level3">
<h3>6.1.1 Number of peaks called</h3>
<pre class="r"><code>##=== R command ===## 
peakN = c()
peakWidth = c()
peakType = c(&quot;control&quot;, &quot;top0.01&quot;)
for(hist in sampleList){
  histInfo = strsplit(hist, &quot;_&quot;)[[1]]
  if(histInfo[1] != &quot;IgG&quot;){
    for(type in peakType){
      peakInfo = read.table(paste0(projPath, &quot;/peakCalling/SEACR/&quot;, hist, &quot;_seacr_&quot;, type, &quot;.peaks.stringent.bed&quot;), header = FALSE, fill = TRUE)  %&gt;% mutate(width = abs(V3-V2))
      peakN = data.frame(peakN = nrow(peakInfo), peakType = type, Histone = histInfo[1], Replicate = histInfo[2]) %&gt;% rbind(peakN, .)
      peakWidth = data.frame(width = peakInfo$width, peakType = type, Histone = histInfo[1], Replicate = histInfo[2])  %&gt;% rbind(peakWidth, .)
    }
  }
}
peakN %&gt;% select(Histone, Replicate, peakType, peakN)</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["chr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["peakType"],"name":[3],"type":["chr"],"align":["left"]},{"label":["peakN"],"name":[4],"type":["int"],"align":["right"]}],"data":[{"1":"K27me3","2":"rep1","3":"control","4":"144906"},{"1":"K27me3","2":"rep1","3":"top0.01","4":"5829"},{"1":"K27me3","2":"rep2","3":"control","4":"74444"},{"1":"K27me3","2":"rep2","3":"top0.01","4":"9642"},{"1":"K4me3","2":"rep1","3":"control","4":"8709"},{"1":"K4me3","2":"rep1","3":"top0.01","4":"878"},{"1":"K4me3","2":"rep2","3":"control","4":"8190"},{"1":"K4me3","2":"rep2","3":"top0.01","4":"3724"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
</div>
<div id="reproducibility-of-the-peak-across-biological-replicates" class="section level3">
<h3>6.1.2 Reproducibility of the peak across biological replicates</h3>
<p>Peak calling on replicate datasets is compared to define reproducible peaks. The top 1% of peaks (ranked by total signal in each block) are selected as high-confidence sites.</p>
<pre class="r"><code>##=== R command ===## 
histL = c(&quot;K27me3&quot;, &quot;K4me3&quot;)
repL = paste0(&quot;rep&quot;, 1:2)
peakType = c(&quot;control&quot;, &quot;top0.01&quot;)
peakOverlap = c()
for(type in peakType){
  for(hist in histL){
    overlap.gr = GRanges()
    for(rep in repL){
      peakInfo = read.table(paste0(projPath, &quot;/peakCalling/SEACR/&quot;, hist, &quot;_&quot;, rep, &quot;_seacr_&quot;, type, &quot;.peaks.stringent.bed&quot;), header = FALSE, fill = TRUE)
      peakInfo.gr = GRanges(peakInfo$V1, IRanges(start = peakInfo$V2, end = peakInfo$V3), strand = &quot;*&quot;)
      if(length(overlap.gr) &gt;0){
        overlap.gr = overlap.gr[findOverlaps(overlap.gr, peakInfo.gr)@from]
      }else{
        overlap.gr = peakInfo.gr
        
      }
    }
    peakOverlap = data.frame(peakReprod = length(overlap.gr), Histone = hist, peakType = type) %&gt;% rbind(peakOverlap, .)
  }
}

peakReprod = left_join(peakN, peakOverlap, by = c(&quot;Histone&quot;, &quot;peakType&quot;)) %&gt;% mutate(peakReprodRate = peakReprod/peakN * 100)
peakReprod %&gt;% select(Histone, Replicate, peakType, peakN, peakReprodNum = peakReprod, peakReprodRate)</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["chr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["peakType"],"name":[3],"type":["chr"],"align":["left"]},{"label":["peakN"],"name":[4],"type":["int"],"align":["right"]},{"label":["peakReprodNum"],"name":[5],"type":["int"],"align":["right"]},{"label":["peakReprodRate"],"name":[6],"type":["dbl"],"align":["right"]}],"data":[{"1":"K27me3","2":"rep1","3":"control","4":"144906","5":"71692","6":"49.47483"},{"1":"K27me3","2":"rep1","3":"top0.01","4":"5829","5":"5314","6":"91.16487"},{"1":"K27me3","2":"rep2","3":"control","4":"74444","5":"71692","6":"96.30326"},{"1":"K27me3","2":"rep2","3":"top0.01","4":"9642","5":"5314","6":"55.11305"},{"1":"K4me3","2":"rep1","3":"control","4":"8709","5":"7959","6":"91.38822"},{"1":"K4me3","2":"rep1","3":"top0.01","4":"878","5":"872","6":"99.31663"},{"1":"K4me3","2":"rep2","3":"control","4":"8190","5":"7959","6":"97.17949"},{"1":"K4me3","2":"rep2","3":"top0.01","4":"3724","5":"872","6":"23.41568"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
<p>The reproducibility is calculated by</p>
<p>`<code># peaks overlapping rep1 and rep2/# peaks of rep1 or rep2 * 100</code></p>
<p>Therefore, it is sensitive to the total number of peaks called in each replicate.</p>
</div>
<div id="fragment-proportion-in-peaks-regions-frips." class="section level3">
<h3>6.1.3 FRagment proportion in Peaks regions (FRiPs).</h3>
<p>We calculate the fraction of reads in peaks (FRiPs) as a measure of signal-to-noise, and contrast it to FRiPs in the IgG control dataset for illustration. Although sequencing depths for CUT&amp;Tag are typically only 1-5 million reads, the low background of the method results in high FRiP scores.</p>
<pre class="r"><code>##=== R command ===## 
library(chromVAR)

bamDir = paste0(projPath, &quot;/alignment/bam&quot;)
inPeakData = c()
## overlap with bam file to get count
for(hist in histL){
  for(rep in repL){
    peakRes = read.table(paste0(projPath, &quot;/peakCalling/SEACR/&quot;, hist, &quot;_&quot;, rep, &quot;_seacr_control.peaks.stringent.bed&quot;), header = FALSE, fill = TRUE)
    peak.gr = GRanges(seqnames = peakRes$V1, IRanges(start = peakRes$V2, end = peakRes$V3), strand = &quot;*&quot;)
    bamFile = paste0(bamDir, &quot;/&quot;, hist, &quot;_&quot;, rep, &quot;_bowtie2.mapped.bam&quot;)
    fragment_counts &lt;- getCounts(bamFile, peak.gr, paired = TRUE, by_rg = FALSE, format = &quot;bam&quot;)
    inPeakN = counts(fragment_counts)[,1] %&gt;% sum
    inPeakData = rbind(inPeakData, data.frame(inPeakN = inPeakN, Histone = hist, Replicate = rep))
  }
}

frip = left_join(inPeakData, alignResult, by = c(&quot;Histone&quot;, &quot;Replicate&quot;)) %&gt;% mutate(frip = inPeakN/MappedFragNum_hg38 * 100)
frip %&gt;% select(Histone, Replicate, SequencingDepth, MappedFragNum_hg38, AlignmentRate_hg38, FragInPeakNum = inPeakN, FRiPs = frip)</code></pre>
<div data-pagedtable="false">
<script data-pagedtable-source type="application/json">
{"columns":[{"label":["Histone"],"name":[1],"type":["chr"],"align":["left"]},{"label":["Replicate"],"name":[2],"type":["chr"],"align":["left"]},{"label":["SequencingDepth"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["MappedFragNum_hg38"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["AlignmentRate_hg38"],"name":[5],"type":["dbl"],"align":["right"]},{"label":["FragInPeakNum"],"name":[6],"type":["dbl"],"align":["right"]},{"label":["FRiPs"],"name":[7],"type":["dbl"],"align":["right"]}],"data":[{"1":"K27me3","2":"rep1","3":"2984630","4":"2859520","5":"95.81","6":"2223813","7":"77.76875"},{"1":"K27me3","2":"rep2","3":"2702260","4":"2606295","5":"96.45","6":"1060373","7":"40.68507"},{"1":"K4me3","2":"rep1","3":"1581710","4":"1494122","5":"94.46","6":"1368641","7":"91.60169"},{"1":"K4me3","2":"rep2","3":"1885056","4":"1742005","5":"92.41","6":"1244857","7":"71.46116"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
</div>
<div id="visualization-of-peak-number-peak-width-peak-reproducibility-and-frips." class="section level3">
<h3>6.1.4 Visualization of peak number, peak width, peak reproducibility and FRiPs.</h3>
<pre class="r"><code>fig7A = peakN %&gt;% ggplot(aes(x = Histone, y = peakN, fill = Histone)) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    facet_grid(~peakType) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.55, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;Number of Peaks&quot;) +
    xlab(&quot;&quot;)

fig7B = peakWidth %&gt;% ggplot(aes(x = Histone, y = width, fill = Histone)) +
    geom_violin() +
    facet_grid(Replicate~peakType) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.55, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    scale_y_continuous(trans = &quot;log&quot;, breaks = c(400, 3000, 22000)) +
    theme_bw(base_size = 18) +
    ylab(&quot;Width of Peaks&quot;) +
    xlab(&quot;&quot;)

fig7C = peakReprod %&gt;% ggplot(aes(x = Histone, y = peakReprodRate, fill = Histone, label = round(peakReprodRate, 2))) +
    geom_bar(stat = &quot;identity&quot;) +
    geom_text(vjust = 0.1) +
    facet_grid(Replicate~peakType) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.55, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;% of Peaks Reproduced&quot;) +
    xlab(&quot;&quot;)

fig7D = frip %&gt;% ggplot(aes(x = Histone, y = frip, fill = Histone, label = round(frip, 2))) +
    geom_boxplot() +
    geom_jitter(aes(color = Replicate), position = position_jitter(0.15)) +
    scale_fill_viridis(discrete = TRUE, begin = 0.1, end = 0.55, option = &quot;magma&quot;, alpha = 0.8) +
    scale_color_viridis(discrete = TRUE, begin = 0.1, end = 0.9) +
    theme_bw(base_size = 18) +
    ylab(&quot;% of Fragments in Peaks&quot;) +
    xlab(&quot;&quot;)

ggarrange(fig7A, fig7B, fig7C, fig7D, ncol = 2, nrow=2, common.legend = TRUE, legend=&quot;bottom&quot;)</code></pre>
<p><img src="figure/tutorials.Rmd/unnamed-chunk-29-1.png" width="100%" style="display: block; margin: auto;" /></p>
<p>
<button type="button" class="btn btn-default btn-xs btn-workflowr btn-workflowr-fig" data-toggle="collapse" data-target="#fig-unnamed-chunk-29-1">
Past versions of unnamed-chunk-29-1.png
</button>
</p>
<div id="fig-unnamed-chunk-29-1" class="collapse">
<div class="table-responsive">
<table class="table table-condensed table-hover">
<thead>
<tr>
<th>
Version
</th>
<th>
Author
</th>
<th>
Date
</th>
</tr>
</thead>
<tbody>
<tr>
<td>
<a href="https://github.com/yezhengSTAT/CUTTag_tutorial/blob/b2993411ac609462d4ebc4b7ce1446e327a00e56/docs/figure/tutorials.Rmd/unnamed-chunk-29-1.png" target="_blank">b299341</a>
</td>
<td>
yezhengSTAT
</td>
<td>
2020-08-10
</td>
</tr>
</tbody>
</table>
</div>
</div>
</div>
</div>
</div>
<div id="vii.-visualization" class="section level1">
<h1>VII. Visualization</h1>
<p>Typically we are interested in visualizing a chromatin landscape in regions using a genome browser. The <a href="http://software.broadinstitute.org/software/igv/home">Integrative Genomic Viewer</a> provides a web app version and a local desktop version that are easy to use. The <a href="https://genome.ucsc.edu/">UCSC Genome Browser</a> provides the most comprehensive supplementary genome information.</p>
<div id="browser-display-of-normalized-bedgraph-files." class="section level2">
<h2>7.1. Browser display of normalized bedgraph files.</h2>
<div class="figure">
<img src="figures/chr7.png" alt="Figure 5. IgV Web Visualization around region chr7:131,000,000-134,000,000" />
<p class="caption"><strong>Figure 5. IgV Web Visualization around region chr7:131,000,000-134,000,000</strong></p>
</div>
</div>
<div id="heatmap-visualization-on-specific-regions" class="section level2">
<h2>7.2. Heatmap visualization on specific regions</h2>
<p>We are also interested in looking at chromatin features at a list of annotated sites, for example histone modification signal at gene promoters.We will use the <code>computeMatrix</code> and <code>plotHeatmap</code> functions from <a href="https://deeptools.readthedocs.io/en/develop/">deepTools</a> to generate the heatmap.</p>
<pre class="bash"><code>##== linux command ==##
mkdir -p $projPath/alignment/bigwig                                                                                                                                        
samtools sort -o $projPath/alignment/bam/${histName}.sorted.bam $projPath/alignment/bam/${histName}_bowtie2.mapped.bam                                                     
samtools index $projPath/alignment/bam/${histName}.sorted.bam                                                                                                              
bamCoverage -b $projPath/alignment/bam/${histName}.sorted.bam -o $projPath/alignment/bigwig/${histName}_raw.bw                                                             
</code></pre>
<div id="heatmap-over-transcription-units" class="section level3">
<h3>7.2.1 Heatmap over transcription units</h3>
<pre class="bash"><code>##== linux command ==##
cores=8
computeMatrix scale-regions -S $projPath/alignment/bigwig/K27me3_rep1_raw.bw \
                               $projPath/alignment/bigwig/K27me3_rep2_raw.bw \
                               $projPath/alignment/bigwig/K4me3_rep1_raw.bw \
                               $projPath/alignment/bigwig/K4me3_rep2_raw.bw \
                              -R $projPath/data/hg38_gene/hg38_gene.tsv \
                              --beforeRegionStartLength 3000 \
                              --regionBodyLength 5000 \
                              --afterRegionStartLength 3000 \
                              --skipZeros -o $projPath/data/hg38_gene/matrix_gene.mat.gz -p $cores

plotHeatmap -m $projPath/data/hg38_gene/matrix_gene.mat.gz -out $projPath/data/hg38_gene/Histone_gene.png --sortUsing sum
</code></pre>
<div class="figure">
<img src="figures/Histone_gene.png" alt="Figure 6. Heatmap of histone enrichment around genes" />
<p class="caption"><strong>Figure 6. Heatmap of histone enrichment around genes</strong></p>
</div>
</div>
<div id="heatmap-on-cuttag-peaks" class="section level3">
<h3>7.2.2. Heatmap on CUT&amp;Tag peaks</h3>
<p>We use the midpoint of the signal block returned from SEACR to align signals in heatmaps. The sixth column of the SEACR output is an entry in the form chr:start-end that represents the first and ending bases of the region with the maximum signal of the region. We first generate a new bed file containing this midpoint information in column 6 and use deeptools for the heatmap visualization.</p>
<pre class="bash"><code>##== linux command ==##
awk &#39;{split($6, summit, &quot;:&quot;); split(summit[2], region, &quot;-&quot;); print summit[1]&quot;\t&quot;region[1]&quot;\t&quot;region[2]}&#39; $projPath/peakCalling/SEACR/${histName}_${repName}_seacr_control.pe\
aks.stringent.bed &gt;$projPath/peakCalling/SEACR/${histName}_${repName}_seacr_control.peaks.summitRegion.bed

computeMatrix reference-point -S $projPath/alignment/bigwig/${histName}_${repName}_raw.bw \
              -R $projPath/peakCalling/SEACR/${histName}_${repName}_seacr_control.peaks.summitRegion.bed \
              --skipZeros -o $projPath/peakCalling/SEACR/${histName}_${repName}_SEACR.mat.gz -p $cores -a 3000 -b 3000 --referencePoint center

plotHeatmap -m $projPath/peakCalling/SEACR/${histName}_SEACR.mat.gz -out $projPath/peakCalling/SEACR/${histName}_SEACR_heatmap.png --sortUsing sum --startLabel &quot;Peak Start&quot; -\
-endLabel &quot;Peak End&quot; --xAxisLabel &quot;&quot; --regionsLabel &quot;Peaks&quot; --samplesLabel &quot;${histName} ${repName}&quot;</code></pre>
<div class="figure">
<img src="figures/Histone_peak.png" alt="Figure 7. Heatmap of histone enrichment in peaks" />
<p class="caption"><strong>Figure 7. Heatmap of histone enrichment in peaks</strong></p>
</div>
</div>
</div>
</div>
<div id="viii.-differential-analysis" class="section level1">
<h1>VIII. Differential analysis</h1>
<ul>
<li>DESeq2: <a href="https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8">Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2</a></li>
</ul>
<p>Estimate variance-mean dependence in count data from high-throughput sequencing assays and test for differential expression based on a model using the negative binomial distribution.</p>
<div id="create-the-peak-x-sample-matrix." class="section level2">
<h2>8.1. Create the peak x sample matrix.</h2>
<p>Usually, the differential tests compare two or more conditions of the same histone modification. In this tutorial, limited by the demonstration data, we will illustrate the differential detection by comparing two replicates of H3K27me3 and two replicates of H3K4me3. We will use DESeq2 (<a href="http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#why-un-normalized-counts">complete tutorial</a>) as illustration.</p>
<div id="create-a-master-peak-list-merging-all-the-peaks-called-for-each-sample." class="section level3">
<h3>8.1.1 Create a master peak list merging all the peaks called for each sample.</h3>
<pre class="r"><code>##=== R command ===## 
mPeak = GRanges()
## overlap with bam file to get count
for(hist in histL){
  for(rep in repL){
    peakRes = read.table(paste0(projPath, &quot;/peakCalling/SEACR/&quot;, hist, &quot;_&quot;, rep, &quot;_seacr_control.peaks.stringent.bed&quot;), header = FALSE, fill = TRUE)
    mPeak = GRanges(seqnames = peakRes$V1, IRanges(start = peakRes$V2, end = peakRes$V3), strand = &quot;*&quot;) %&gt;% append(mPeak, .)
  }
}
masterPeak = reduce(mPeak)</code></pre>
</div>
<div id="get-the-fragment-counts-for-each-peak-in-the-master-peak-list." class="section level3">
<h3>8.1.2 Get the fragment counts for each peak in the master peak list.</h3>
<pre class="r"><code>##=== R command ===## 
library(DESeq2)
bamDir = paste0(projPath, &quot;/alignment/bam&quot;)
countMat = matrix(NA, length(masterPeak), length(histL)*length(repL))
## overlap with bam file to get count
i = 1
for(hist in histL){
  for(rep in repL){
    
    bamFile = paste0(bamDir, &quot;/&quot;, hist, &quot;_&quot;, rep, &quot;_bowtie2.mapped.bam&quot;)
    fragment_counts &lt;- getCounts(bamFile, masterPeak, paired = TRUE, by_rg = FALSE, format = &quot;bam&quot;)
    countMat[, i] = counts(fragment_counts)[,1]
    i = i + 1
  }
}
colnames(countMat) = paste(rep(histL, 2), rep(repL, each = 2), sep = &quot;_&quot;)</code></pre>
</div>
</div>
<div id="sequencing-depth-normalization-and-differential-enriched-peaks-detection" class="section level2">
<h2>8.2. Sequencing depth normalization and differential enriched peaks detection</h2>
<pre class="r"><code>##=== R command ===## 
selectR = which(rowSums(countMat) &gt; 5) ## remove low count genes
dataS = countMat[selectR,]
condition = factor(rep(histL, each = length(repL)))
dds = DESeqDataSetFromMatrix(countData = dataS,
                              colData = DataFrame(condition),
                              design = ~ condition)
DDS = DESeq(dds)
normDDS = counts(DDS, normalized = TRUE) ## normalization with respect to the sequencing depth
colnames(normDDS) = paste0(colnames(normDDS), &quot;_norm&quot;)
res = results(DDS, independentFiltering = FALSE, altHypothesis = &quot;greaterAbs&quot;)

countMatDiff = cbind(dataS, normDDS, res)
head(countMatDiff)</code></pre>
<pre><code>DataFrame with 6 rows and 14 columns
  K27me3_rep1 K4me3_rep1 K27me3_rep2 K4me3_rep2 K27me3_rep1_norm
    &lt;numeric&gt;  &lt;numeric&gt;   &lt;numeric&gt;  &lt;numeric&gt;        &lt;numeric&gt;
1           6          2           1          6         1.408657
2           1          0         242        182         0.234776
3           0          0         176         88         0.000000
4           0          0         274        194         0.000000
5           3          4           0          1         0.704328
6           0          1         109         59         0.000000
  K4me3_rep1_norm K27me3_rep2_norm K4me3_rep2_norm  baseMean log2FoldChange
        &lt;numeric&gt;        &lt;numeric&gt;       &lt;numeric&gt; &lt;numeric&gt;      &lt;numeric&gt;
1        0.620403            4.170         18.2724   6.11787       3.496854
2        0.000000         1009.141        554.2634 390.90978      12.510325
3        0.000000          733.921        267.9955 250.47905      13.297304
4        0.000000         1142.581        590.8082 433.34733      14.089840
5        1.240806            0.000          3.0454   1.24763       0.846266
6        0.310202          454.530        179.6788 158.62986      11.189689
      lfcSE      stat      pvalue        padj
  &lt;numeric&gt; &lt;numeric&gt;   &lt;numeric&gt;   &lt;numeric&gt;
1   1.19893  2.916635 3.53829e-03 4.22134e-02
2   1.50039  8.338074 7.55102e-17 2.18197e-15
3   1.58547  8.386969 4.98837e-17 1.50548e-15
4   1.55196  9.078730 1.09850e-19 6.73560e-18
5   2.18326  0.387617 6.98300e-01 9.72755e-01
6   1.53046  7.311313 2.64545e-13 4.33546e-12</code></pre>
<ul>
<li><p>DESeq2 requires the input matrix should be un-normalized counts or estimated counts of sequencing reads.</p></li>
<li><p>DESeq2 model internally corrects for library size.</p></li>
<li><code>countMatDiff</code> summarizes the differential analysis results:
<ul>
<li>First 4 columns: raw reads counts after filtering the peak regions with low counts</li>
<li>Second 4 columns: normalized read counts eliminating library size difference.</li>
<li>Remaining columns: differential detection results.</li>
</ul></li>
</ul>
</div>
</div>
<div id="references" class="section level1">
<h1>References</h1>
<pre><code>Kaya-Okur HS, Wu SJ, Codomo CA, Pledger ES, Bryson TD, Henikoff JG, Ahmad K, Henikoff S: CUT&amp;Tag for efficient epigenomic profiling of small samples and single cells. Nature Communications 2019 10:1930 (PMID:31036827).</code></pre>
<pre><code>Meers, M.P., Tenenbaum, D. &amp; Henikoff, S. Peak calling by Sparse Enrichment Analysis for CUT&amp;RUN chromatin profiling. Epigenetics &amp; Chromatin 12, 42 (2019). https://doi.org/10.1186/s13072-019-0287-4</code></pre>
</div>
<div id="cite-this-tutorial" class="section level1">
<h1>Cite this tutorial</h1>
<p>Zheng Y et al (2020). Protocol.io</p>
</div>
<div id="ix.-additional-alternatives" class="section level1">
<h1>IX. Additional Alternatives</h1>
<div id="chipseqspikeinfree-for-normalizing-data-without-spike-in-dna-optional" class="section level2">
<h2>9.1 ChIPseqSpikeInFree for normalizing data without spike-in DNA [Optional]</h2>
<p><a href="https://academic.oup.com/bioinformatics/article/36/4/1270/5578481">ChIPseqSpikeInFree: a ChIP-seq normalization approach to reveal global changes in histone modifications without spike-in</a> is a novel ChIP-seq normalization method to effectively determine scaling factors for samples across various conditions and treatments, which does not rely on exogenous spike-in chromatin or peak detection to reveal global changes in histone modification occupancy. The installation details can be found on <a href="https://github.com/stjude/ChIPseqSpikeInFree">github</a>.</p>
<ul>
<li><p><a href="https://github.com/stjude/ChIPseqSpikeInFree#interpretation-of-scaling-factor-table">Interpretation of the ChIPseqSpikeInFree output.</a></p></li>
<li><p><a href="https://github.com/stjude/ChIPseqSpikeInFree#how-to-use-chipseqspikein-scaling-factor">How to use ChIPseqSpikeInFree scaling factor.</a></p></li>
</ul>
</div>
<div id="other-peak-calling-methods." class="section level2">
<h2>9.2. Other peak calling methods.</h2>
<ul>
<li>MACS2: <a href="https://genomebiology.biomedcentral.com/articles/10.1186/gb-2008-9-9-r137">Model-based Analysis of ChIP-Seq (MACS)</a>. Installation details can be found <a href="https://github.com/taoliu/MACS/wiki">here</a>.</li>
</ul>
<pre class="bash"><code>##== linux command ==##
histName=&quot;K27me3&quot;
controlName=&quot;IgG&quot;

mkdir -p $projPath/peakCalling
macs2 callpeak -t ${projPath}/alignment/bam/${histName}_rep1_bowtie2.mapped.bam \
      -c ${projPath}/alignment/bam/${controlName}_rep1_bowtie2.mapped.bam \
      -g hs -f BAMPE -n macs2_peak_q0.1 --outdir $projPath/peakCalling/MACS2 -q 0.1 --keep-dup all 2&gt;${projPath}/peakCalling/MACS2/macs2Peak_summary.txt</code></pre>
<ul>
<li><p>dPeak: <a href="https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003246">dPeak: High Resolution Identification of Transcription Factor Binding Sites from PET and SET ChIP-Seq Data</a></p></li>
<li><p>MOSAiCS: <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4608541/">A Statistical Framework for the Analysis of ChIP-Seq Data</a></p></li>
</ul>
</div>
<div id="other-packages-for-differential-analysis-of-binding-sites" class="section level2">
<h2>9.3 Other packages for differential analysis of binding sites</h2>
<ul>
<li>Limma: <a href="https://academic.oup.com/nar/article/43/7/e47/2414268">limma powers differential expression analyses for RNA-sequencing and microarray studies</a></li>
</ul>
<p>Limma is an R package for the analysis of gene expression microarray data, especially the use of linear models for analysing designed experiments and the assessment of differential expression. Limma provides the ability to analyse comparisons between many RNA targets simultaneously in arbitrary complicated designed experiments. Empirical Bayesian methods are used to provide stable results even when the number of arrays is small. Limma can be extended to study differential fragment enrichment analysis within peak regions. Notably, limma can deal with both fixed effect model and random effect model.</p>
<ul>
<li>edgeR: <a href="https://academic.oup.com/nar/article/40/10/4288/2411520">Differential Expression Analysis of Multifactor RNA-Seq Experiments With Respect to Biological Variation</a></li>
</ul>
<p>Differential expression analysis of RNA-seq expression profiles with biological replication. Implements a range of statistical methodology based on the negative binomial distributions, including empirical Bayes estimation, exact tests, generalized linear models and quasi-likelihood tests. As well as RNA-seq, it be applied to differential signal analysis of other types of genomic data that produce read counts, including ChIP-seq, ATAC-seq, Bisulfite-seq, SAGE and CAGE. edgeR can deal with multifactor problem.</p>
</div>
</div>
<div id="x.-troubleshooting-generating-your-data" class="section level1">
<h1>X. Troubleshooting: Generating your data</h1>
<p>This workflow can be followed with your own data and will generate a standardized set of quality-control reports. However, many sequencing facilities do not perform 25x25 PE sequencing, and alternate parameters for trimming and mapping are provided here. Control datasets for non-specific antibody (IgG) profiling or ATAC-seq profiling of your material can also be used for optional analysis detailed here.</p>
<p>Stringent washing with 300 mM NaCl is critical to limit the affinity of Tn5 for exposed DNA. We describe here the need for controlling background Tn5 affinities and describe how our CUT&amp;Tag protocol effectively suppresses this artifact for unambiguous mapping of chromatin epitopes. We present a protocol that can process either native or fixed nuclei and includes alternative methods for DNA isolation. To illustrate the method, we describe a typical experiment, including evaluation of the results using a new metric for peak-calling information. Further, we validate a single-tube format for CUT&amp;Tag that requires no DNA isolation but instead uses tagmented material directly for library amplification. We document critical steps for the CUT&amp;Tag protocol, informed by our experiences, helping users establish this method in their research.</p>
<br>
<p>
<button type="button" class="btn btn-default btn-workflowr btn-workflowr-sessioninfo" data-toggle="collapse" data-target="#workflowr-sessioninfo" style="display: block;">
<span class="glyphicon glyphicon-wrench" aria-hidden="true"></span> Session information
</button>
</p>
<div id="workflowr-sessioninfo" class="collapse">
<pre class="r"><code>sessionInfo()</code></pre>
<pre><code>R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS/LAPACK: /app/software/OpenBLAS/0.3.7-GCC-8.3.0/lib/libopenblas_haswellp-r0.3.7.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] corrplot_0.84               ggpubr_0.4.0               
 [3] DESeq2_1.28.1               SummarizedExperiment_1.18.2
 [5] DelayedArray_0.14.1         matrixStats_0.57.0         
 [7] Biobase_2.48.0              chromVAR_1.10.0            
 [9] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2        
[11] IRanges_2.22.2              S4Vectors_0.26.1           
[13] BiocGenerics_0.34.0         viridis_0.5.1              
[15] viridisLite_0.3.0           ggplot2_3.3.2              
[17] stringr_1.4.0               dplyr_1.0.2                
[19] workflowr_1.6.2            

loaded via a namespace (and not attached):
  [1] colorspace_1.4-1            ggsignif_0.6.0             
  [3] rio_0.5.16                  ellipsis_0.3.1             
  [5] rprojroot_2.0.2             XVector_0.28.0             
  [7] fs_1.4.2                    rstudioapi_0.11            
  [9] farver_2.0.3                DT_0.16                    
 [11] bit64_4.0.5                 AnnotationDbi_1.50.3       
 [13] splines_4.0.2               R.methodsS3_1.8.0          
 [15] geneplotter_1.66.0          knitr_1.29                 
 [17] jsonlite_1.7.2              Cairo_1.5-12               
 [19] Rsamtools_2.4.0             seqLogo_1.54.3             
 [21] broom_0.5.6                 annotate_1.66.0            
 [23] GO.db_3.11.4                png_0.1-7                  
 [25] R.oo_1.23.0                 shiny_1.5.0                
 [27] readr_1.3.1                 compiler_4.0.2             
 [29] httr_1.4.2                  backports_1.2.1            
 [31] Matrix_1.2-18               fastmap_1.0.1              
 [33] lazyeval_0.2.2              later_1.1.0.1              
 [35] htmltools_0.5.0             tools_4.0.2                
 [37] gtable_0.3.0                glue_1.4.2                 
 [39] TFMPvalue_0.0.8             GenomeInfoDbData_1.2.3     
 [41] reshape2_1.4.4              Rcpp_1.0.5                 
 [43] carData_3.0-4               cellranger_1.1.0           
 [45] vctrs_0.3.5                 Biostrings_2.56.0          
 [47] nlme_3.1-148                rtracklayer_1.48.0         
 [49] xfun_0.15                   CNEr_1.24.0                
 [51] openxlsx_4.1.5              mime_0.9                   
 [53] miniUI_0.1.1.1              lifecycle_0.2.0            
 [55] poweRlaw_0.70.6             gtools_3.8.2               
 [57] rstatix_0.6.0               XML_3.99-0.5               
 [59] zlibbioc_1.34.0             scales_1.1.1               
 [61] BSgenome_1.56.0             hms_0.5.3                  
 [63] promises_1.1.1              RColorBrewer_1.1-2         
 [65] curl_4.3                    yaml_2.2.1                 
 [67] memoise_1.1.0               gridExtra_2.3              
 [69] stringi_1.5.3               RSQLite_2.2.0              
 [71] highr_0.8                   genefilter_1.70.0          
 [73] caTools_1.18.0              zip_2.0.4                  
 [75] BiocParallel_1.22.0         rlang_0.4.9                
 [77] pkgconfig_2.0.3             bitops_1.0-6               
 [79] pracma_2.2.9                evaluate_0.14              
 [81] lattice_0.20-41             purrr_0.3.4                
 [83] labeling_0.4.2              GenomicAlignments_1.24.0   
 [85] htmlwidgets_1.5.3           cowplot_1.1.0              
 [87] bit_4.0.4                   tidyselect_1.1.0           
 [89] plyr_1.8.6                  magrittr_1.5               
 [91] R6_2.5.0                    generics_0.0.2             
 [93] DBI_1.1.0                   haven_2.3.1                
 [95] foreign_0.8-80              pillar_1.4.7               
 [97] whisker_0.4                 withr_2.3.0                
 [99] abind_1.4-5                 survival_3.2-3             
[101] KEGGREST_1.28.0             RCurl_1.98-1.2             
[103] tibble_3.0.4                car_3.0-8                  
[105] crayon_1.3.4                plotly_4.9.2.1             
[107] rmarkdown_2.3               TFBSTools_1.26.0           
[109] readxl_1.3.1                locfit_1.5-9.4             
[111] grid_4.0.2                  data.table_1.13.4          
[113] blob_1.2.1                  git2r_0.27.1               
[115] forcats_0.5.0               digest_0.6.27              
[117] xtable_1.8-4                tidyr_1.1.2                
[119] httpuv_1.5.4                R.utils_2.9.2              
[121] munsell_0.5.0               DirichletMultinomial_1.30.0</code></pre>
</div>
</div>


<!-- Adjust MathJax settings so that all math formulae are shown using
TeX fonts only; see
http://docs.mathjax.org/en/latest/configuration.html.  This will make
the presentation more consistent at the cost of the webpage sometimes
taking slightly longer to load. Note that this only works because the
footer is added to webpages before the MathJax javascript. -->
<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    "HTML-CSS": { availableFonts: ["TeX"] }
  });
</script>


</div>
</div>

</div>

<script>

// add bootstrap table styles to pandoc tables
function bootstrapStylePandocTables() {
  $('tr.header').parent('thead').parent('table').addClass('table table-condensed');
}
$(document).ready(function () {
  bootstrapStylePandocTables();
});


</script>

<!-- tabsets -->

<script>
$(document).ready(function () {
  window.buildTabsets("TOC");
});

$(document).ready(function () {
  $('.tabset-dropdown > .nav-tabs > li').click(function () {
    $(this).parent().toggleClass('nav-tabs-open')
  });
});
</script>

<!-- code folding -->
<script>
$(document).ready(function () {
  window.initializeCodeFolding("show" === "show");
});
</script>

<script>
$(document).ready(function ()  {

    // move toc-ignore selectors from section div to header
    $('div.section.toc-ignore')
        .removeClass('toc-ignore')
        .children('h1,h2,h3,h4,h5').addClass('toc-ignore');

    // establish options
    var options = {
      selectors: "h1,h2,h3",
      theme: "bootstrap3",
      context: '.toc-content',
      hashGenerator: function (text) {
        return text.replace(/[.\\/?&!#<>]/g, '').replace(/\s/g, '_');
      },
      ignoreSelector: ".toc-ignore",
      scrollTo: 0
    };
    options.showAndHide = true;
    options.smoothScroll = true;

    // tocify
    var toc = $("#TOC").tocify(options).data("toc-tocify");
});
</script>

<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
